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We perform a new analysis of electron-proton scattering data to determine the proton electric 
and magnetic radii, enforcing model-independent constraints from form factor analyticity. A wide- 
ranging study of possible systematic effects is performed. An improved analysis is developed that 
rebins data taken at identical kinematic settings, and avoids a scaling assumption of systematic 

errors with statistical errors. Employing standard models for radiative corrections, our improved 

analysis of the 2010 Mainz A1 Collaboration data yields a proton electric radius te = 0.895(20) fm 
and magnetic radius vm = 0.776(38) fm. A similar analysis applied to world data (excluding 
Mainz data) implies rE = 0.916(24) fm and vm = 0.914(35) fm. The Mainz and world values 

of the charge radius are consistent, and a simple combination yields a value r\e = 0.904(15) fm 

that is 4 <j larger than the CREMA Collaboration muonic hydrogen determination. The Mainz and 
world values of the magnetic radius differ by 2.7 ct, and a simple average yields tm = 0.851(26) fm. 
The circumstances under which published muonic hydrogen and electron scattering data could be 
reconciled are discussed, including a possible deficiency in the standard radiative correction model 
which requires further analysis. 
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I. INTRODUCTION 


The electromagnetic form factors of the nucleons pro¬ 
vide basic inputs to precision tests of the Standard Model 
and to the determination of fundamental constants pQ. 
These form factors are also of critical importance for the 
accelerator neutrino program [2]. The development of 
muonic atom spectroscopy has introduced a power¬ 
ful new probe of proton and nuclear structure, challeng¬ 
ing existing results obtained from (electronic) hydrogen 
and electron scattering Q] . Taken at face value, in the 
absence of new physics explanations, the muonic hydro¬ 
gen Lamb shift measurement [3] necessitates a > 5cr re¬ 
vision of the Rydberg constant, in addition to discarding 
or revising the predictions from a large body of previous 
results in both electron-proton scattering and hydrogen 
spectroscopy. Sources of systematic error in electron- 
proton scattering measurements also impact neutrino- 
nucleus scattering observables and hence the extraction 
of fundamental neutrino oscillation parameters at cur¬ 
rent and future facilities mw- Resolution of the so- 
called proton radius puzzle thus has important implica¬ 
tions across the fields of high energy, nuclear, and atomic 
physics [TUB]. 

The muonic hydrogen measurement [J yields r# = 
0.84087(39) fm, compared to t\e = 0.8758(77) fm for 
Lamb shift measurements from ordinary (electronic) hy¬ 
drogen [lj. Previous analyses of electron scattering re¬ 
sults using the high statistics data set taken at the Mainz 
Microtron (MAMI) yielded HiHTT] te = 0.879(11) fm 
and vm = 0.777(19) fm, in both cases neglecting un¬ 
certainty associated with two-photon exchange correc¬ 
tions m- A previous global analysis of world data m , 
excluding the Mainz data, yielded t\e = 0.875(10) fm 
and vm = 0.867(20) fm. Similar results were obtained 
in an independent global analysis which included con¬ 
straints on the large-distance behavior of an inferred pro¬ 
ton charge distribution mm- So not only are electron- 
and muon-based extractions of the charge radius incon¬ 
sistent, but there is also a ~ 3<r disagreement between 
extractions of tm from different electron scattering data 
sets. 

Here, we address the issue of radius extraction from 
electron-proton scattering data. A prominent uncer¬ 
tainty in the extracted radius arises from the shape of 
the form factor assumed when extrapolating to q 2 = 0 
where the radius is defined in terms of the form factor 
slope. This can be the dominant uncertainty, as hap¬ 
pens in particular in the A1 Collaboration’s extraction of 
the charge radius from Mainz data [S]. In Ref. m , one 
of us investigated the implications of analyticity for the 
form factors of electromagnetic lepton-nucleon scatter¬ 
ing. Reference considered a representative data set 
consisting of extracted electric form factor values from 
cross section data prior to 2007. In the present paper, 
we extend this analysis by fitting directly to cross sec¬ 
tion data, eliminating possible systematic uncertainties 
associated with the reduction from cross sections to form 


factors prior to the q 2 —» 0 extrapolation that defines the 
radius observable. We consider the most recent data, in¬ 
cluding separately a “Mainz dat a set” [3] and a “world 
data set” (defined below in Sec. IV E| excluding Mainz 
data). We extend the analysis of the electric form factor 
of the proton to consider also the magnetic form factor 
(see also Ref. [17]), necessary to connect with cross sec¬ 
tion data. We discuss the impact of uncertainties arising 
from the fitting procedure, theoretical corrections to the 
cross sections, and experimental systematic uncertain¬ 
ties. We focus here on understanding the implications 
of electron-proton scattering data in isolation and do not 
include further constraints arising from isospin decompo¬ 
sition in combination with electron-neutron, 7T7t —> NN, 
or other data [TTH - lTSl . As we will see, several critical 
issues in the electron scattering data demand attention 
before the inclusion of such further ingredients. 

The remainder of the paper is structured as follows. 
Section mi introduces notations and conventions. Sec¬ 
tion III discusses the constraints of the form factor shape 
arising from analyticity and perturbative scaling. Sec¬ 
tion m reviews the status of radiative corrections and 
defines the default models used in the remainder of the 
paper. Section |V| analyzes the Mainz data set, employing 
exactly the same analysis strategy as detailed in Ref. [9], 
with the exception that the bounded z expansion is used 
in place of polynomial or spline functions to represent 
the form factors. Section VI studies a range of possible 
systematic effects. Section VII| provides updated extrac¬ 
tions of the charge and magnetic radii and uncertainties, 
for both the Mainz data and the world data. Section lVlHI 
presents a summary and conclusions. Supplemental Ma¬ 
terial [13 includes the data used for fits in Sections [V VI 
and IVIII 


II. CONVENTIONS AND NOTATION 

The Dirac and Pauli form factors of the proton, i 7 ) and 
F‘ 2 , respectively, are defined by 

{p(p')\ J£n\p(p)) = u(p')r (p) »(p',p)u(p) , (1) 

where 

r (p) ») = 7^1 (9 2 ) + 2^ q 2 ), (2) 

with q p = p'^—p^. The Sachs electric and magnetic form 
factors are related to the Dirac-Pauli basis by 

G E (q 2 ) = F 1 (q 2 ) + ^F 2 (q 2 ), 

G M (q 2 ) = F^q 2 ) + F 2 (q 2 ), (3) 

where Ge( 0) = 1, Gm{ 0) = p p « 2.793 J2U]. The electric 
and magnetic radii, re and rjf, are defined as the slopes 
of the Sachs form factors at q 2 = 0, i.e., 

G_b,m(9 2 ) , 1 2 2 , 4\ fA\ 

GeM 0 ) " G E ’ Mq ( ’ q) ' ( } 
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In terms of G E and Gm, the cross section for electron- 
proton scattering in single photon exchange approxima¬ 
tion is 

fda\ = fda\ eG 2 E +rG 2 M 

\dn ) 0 \dnj Mott e ( 1 + ' r ) 

where (da/dfl)Mott is the recoil-corrected relativistic 
point particle (Mott) result, 


/ da \ a 2 E' 2 9 

\dn) M ott = 4^W|^ C0S '2- 


( 6 ) 


Here, Q 2 = —q 2 , E is the initial electron energy, E' = 
E/[l + (2 E/rrip) sin 2 (#/2)] and 9 are the energy and an¬ 
gle with respect to the beam direction of the final state 
electron, and e, r are the dimensionless kinematic vari¬ 
ables 


Q 2 

4 m 2 ’ 


e 


1 + 2(1 + r) tan 2 


9' 

2 



(7) 


In fits to the Mainz data below, we employ the beam 
energy E and the acceptance-averaged Q 2 as independent 
variables, as dictated by the presentation of experimental 
results in this data set. In fits to world data excluding 
Mainz data, we employ E and 9. 


a given experimental data set, — Qmax < t < 0, the finite 
distance to singularities implies the existence of a small 
expansion parameter, |z| max < 1. To see this, we perform 
a conformal mapping of the domain of analyticity onto 
the unit circle, 2 


z(t,t cut ,t 0 ) = 


V^cut t \/ tent. to 
\J G:ut t -f- y^cut to 


( 8 ) 


where t C ut = 4m 2 and to is a free parameter repre¬ 
senting the point mapping onto 2 = 0. By the “op¬ 
timal” choice to Pt (Qmax) = *cut (l - \/l + QmaxAcut) , 
the maximum value of \z\ is minimized: |^| < l^lSfax = 

[(! + QLxAcut)* - 1]/[(1 + <5LxAcut)^ + 1]■ For exam¬ 
ple, with Q 2 lax = 0.05, 0.5, and 1 GeV 2 , we have |2|^ x ~ 
0.06, 0.25, and 0.32, respectively. 3 Expanding the form 
factors as 


G E (q 2 )= G M (q 2 )=Y, bkzk > ( 0 ) 

k =0 fc=0 


we find that higher-order terms are suppressed by powers 
of this small parameter. 


B. Coefficient bounds and large-fc scaling 


III. FORM FACTOR SHAPE 

When performing statistical analyses that constrain 
the form factors and derived quantities such as the ra¬ 
dius, it is important that the class of allowed functions 
be large enough to contain the true form factors but suffi¬ 
ciently constrained for meaningful values to be obtained, 
i.e., without arbitrarily large errors, and such that over¬ 
fitting to statistical noise does not bias parameter ex¬ 
tractions. We summarize here our knowledge about the 
analytic structure of the form factors, introduce notation 
for the z expansion, and explain an important property 
of the 2 expansion with regard to convexity and stability 
of fits involving many parameters. 


A. Analyticity and 2 expansion 

Let us recall the analytic structure of the form factors 
F 1 (q 2 ), F 2 (q 2 ), or equivalently G E {q 2 ), G M (q 2 )- The 
form factors may be extended to functions of the complex 
variable t = q 2 , analytic outside of a cut at timelike val¬ 
ues of t, beginning at the two-pion production threshold, 
t > 4 m 2 . 1 In the restricted kinematic region accessed in 


The identity |1|, 



k =o 


1 dt 

7T 7 tout t - t 0 



\G\ 2 <oo, 


( 10 ) 


ensures that the coefficients multiplying z k are not only 
bounded in size but must decrease at large k. This guar¬ 
antees that a finite number of parameters is necessary to 
describe the form factor with a given precision through¬ 
out the kinematic region of interest. In the fits per¬ 
formed in this paper, we focus on the class of form fac¬ 
tors with a uniform bound on \ak/ao\ and \bk/bo\. A 
study of form factor models, explicit spectral functions, 
and scattering data motivates the conservative bound of 
|a fc /a 0 |max = |&fc/&o|max = 5 for either t 0 = 0 or t 0 = t^ pt 
when limiting the fit to Q^ ax « 1 GeV 2 mm We adopt 
this bound for our default fits but study also the case of 
modified bounds. 

In fact, a stronger statement can be made regarding 
the large-A; behavior of the expansion coefficients. Since 
at large spacelike values of momentum transfer, Q 2 —> 
oo, the Sachs form factors are known to fall as 1 /Q 4 
up to logarithms (223, we have that Q k G(—Q 2 ) —> 0, 
k = 0,..., 3. From Eq. (|9j) this implies 


d n 

dz n 


G e 


= 0, 

Z=1 


n = 0,1, 2, 3, 


( 11 ) 


1 Here and throughout, = 140 MeV denotes the charged pion 
mass. Accounting for isospin violation would imply a smaller 
threshold at 4m^ 0 . A conservative approach to accounting for 
this effect would be to lower the threshold; we have verified that 
the difference is inconsequential to the fits. 


2 For a discussion and further references, see Ref. m • 

3 For to — 0, these numbers become approximately twice as large, 

| z |max ~ 0.12, 0.46, and 0.58. 
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or, equivalently, the series of four sum rules, 

OO 

k(k — 1) • • • (k — n + 1 )afc = 0, n = 0,1,2,3. (12) 

k—n 

Absolute convergence of these series corresponds to a & = 
C>(l/fc 4 ). 4 


C. Convexity and y 2 minimization 


The preceding arguments strictly apply when a single 
form factor dominates, which is, for example, the case 
for Ge in Eq. ([ 5 ]) at low Q 2 . In the general case in which 
Ge and Gm are fit simultaneously, the y 2 function takes 
a more general form involving the sum of probabilities, 
A 2 —> A 2 + Bf, for which the simple convexity theorem 
following from Eq. (14) no longer applies. It may be 
interesting to pursue more general “physical convexity” 
theorems involving multiple probability sums and corre¬ 
lated errors. 


An important feature of Taylor expanded amplitudes 
facilitates efficient and stable numerical fits. 5 Consider a 
X 2 function of schematic form 

X 2 = E (A1 ‘ )2 , A = £«"*."• (13) 

i 1 n 


where the sum is over data points labeled by index i, M 
represents a measurement, E is the error on the mea¬ 
surement, and A is the theoretical amplitude expressed 
in terms of a kinematic variable x [e.g., Xi = z(q 2 ) in the 
present application]. The Taylor expansion coefficients 
a n are fit parameters to be determined by minimizing 

x 2 - 

If X 2 ({°n}) is a convex function of its arguments, a„, 
then any local minimum is necessarily a global mini¬ 
mum, and the relevant optimization problem is amenable 
to efficient numerical algorithms. In general, determin¬ 
ing convexity of a multivariate quartic polynomial is NP 
hard [22]. We notice, however, that the matrix of second 
derivatives is 


d 2 X 2 = n +m ( 3 A 2 - M.j \ 

da n da m ^ 1 \ Ef ) ■ 


(14) 


Each term in the sum over i is seen to be a posi¬ 
tive semidefinite matrix provided that a mpl itudes satisfy 
3 A 2 > Mi. Each contribution in Eq. (13) is thus con¬ 
vex throughout the parameter regime where this physical 
condition is satisfied. Since a linear combination of con¬ 
vex polynomials with positive coefficients is convex, the 
sum of terms in Eq. (13) is also convex in this regime. 
It is thus straightforward to build up solutions to the 
numerical y 2 minimization problem over a large number 
of parameters by successively increasing the number of 
parameters and data points, using the previous solution 
{a n } as the initial condition. The convexity condition en¬ 
sures that this procedure does not yield a solution that 
is a local but not a global minimum. 


4 Absolu te c onvergence may be verified by inspecting the analog 
of Eq. jlol applied to \d n G/dz n \ 2 in place of |G| 2 . 

5 While this observation has emerged from a particular example 

of fits to electron scattering data, the argument applies to gen¬ 
eral quantum mechanical observables represented as squares of 
Taylor-expanded amplitudes. 


D. Advantages over other parametrizations 


We remark that several parametrizations of the pro¬ 
ton form factors in common use rely on somewhat arbi¬ 
trary expansions. A simple Taylor expansion in q 2 [24] 
is only guaranteed to converge below the pion produc¬ 
tion threshold q 2 < 4 m 2 « 0.08 GeV 2 . Convergence of a 
sequence of Pade approximants, implemented either di¬ 
rectly as a ratio of polynomials [25], (2B] or as a continued 
fraction [27] [28], requires positivity of the spectral func¬ 
tion in the dispersive representations of the form factors, 
a property which is not satisfied. 6 

While these functions may be able to provide a suf¬ 
ficiently precise representation of the form factors with 
enough fit parameters, the parameters tend to be highly 
correlated. Without any way to bound the parameters, 
these correlations can lead to a large uncertainty on any 
given parameter (such as the radius) that grows as the 
number of parameters increases. Because of this, it may 
be difficult or impossible to include enough parameters 
to properly reproduce the data while at the same time 
achieving a meaningful limit on the extracted radius. The 
correlation between different parameters may also lead to 
the situation in which overfitting the noise in data at high 
Q 2 biases the extracted radius. This concern applies es¬ 
pecially for the magnetic form factor for which the data 
at low Q 2 have larger uncertainties than the higher Q 2 
data. 


IV. RADIATIVE CORRECTIONS 


In this section, we provide a brief summary of one- 
loop radiative corrections and Sudakov resummation in 
electron-proton scattering. We extract the radius accord¬ 
ing to Eq. Q , using data to which corrections have been 
applied to extract a Born cross section. To understand 
the impact of the different corrections applied to vari¬ 
ous data sets, we begin in Secs. IV A| |IV D[ with a brief 
overview of notation and results for one-photon exchange, 


That it cannot be satisfied is readily seen from the asymptotic 
behavior Q~ 2 for the form factor represented by such a spectral 
function. 
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two-photon exchange, real photon emission, and Sudakov 
resummation as they impact cross section measurements. 
In Sec. |IV E[ we return to a discussion of experimental 
implementations for the data sets employed in the re¬ 
mainder of this paper. 


A. One-photon exchange 


The (on-shell, renormalized) scattering amplitude 
for the electron-proton scattering process e(k)p(p) —> 
e{k')p{p') involving one exchanged photon may be writ¬ 
ten 


47TQf 


1 


M 1 = t u^\k')T^{k',k)u^{k) 

r l-n (q 2 ) 

x u (p) (p , ) r ^ ) (p',p)u {p) (p) , (15) 

and includes radiative corrections involving the proton 
and electron vertices (and wave function renormaliza¬ 
tion) and vacuum polarization. Here, a = 7.297 x 10“ 3 
is the fine structure constant. The proton vertex func¬ 
tion T^ p \p',p) is expressed, as in Eq. in terms of 
the IR divergent on-shell form factors discussed below. 
The electron vertex function T^ e \k\k) is similarly ex¬ 
pressed in terms of on-shell form factors normalized as 
F^ e \ 0) = 1, F^fi) = a e ~ a/(2w). The photon propa¬ 
gator correction II(< 7 2 ) accounts for contributions of both 
leptonic and hadronic vacuum polarization. 


The on-shell form factors appearing in Eq. (151 are nec¬ 


essarily infrared divergent at nonzero momentum trans¬ 
fer, as deduced by the cancellation with bremsstrahlung 
emission. In terms of a photon mass, let us introduce 
conventional “Born” form factors which are finite includ¬ 
ing first-order radiative effects in the A —> 0 limit. We 
employ the tilde notation Fj to denote the on-shell form 
factor with the corresponding Born form factor F,: 

Fi(<i 2 ) = j 1 - ^ [K{p,p') - K(p,p )] | Ffiq 2 ). (16) 

Here K(p 1 ,p 2 ) denotes the integral [25] {301 

1 


Tf, \ 2pi -p 2 /■ 4 

K (Pi , P2) = -/ d L 


— Z7T 


1 


L 2 - A 2 + iO 
1 


L 2 + 2L ■ pi + iO L 2 + 2L ■ p 2 + iO 


(17) 


and is readily evaluated in analytic form. The electric 
and magnetic radii are now defined as the slopes of the 
Born form factors with respect to q 2 , 




G e ,m( 0) 6 rE ' M ' 


(18) 


Infrared divergences are absorbed into the extracted pref¬ 
actors in Eq. (16) and will cancel upon including the ef¬ 


fects of real photon emission. The electron form factors 


may be calculated analytically in QED, with infrared di¬ 
vergences similarly cancelling against real emission. For 
completeness, let us note that for the IR divergent on- 
shell Sachs form factors we have 


_ 1^2 

Ge,m{ 0 ) 6 


a 

37t m 2 




(19) 


Leptonic vacuum polarization contributions to n(g 2 ) 
are readily computed analytically, and the hadronic con¬ 
tributions are constrained by e + e~ —> hadrons data. 7 
For the purposes of determining the radii, we may sim¬ 
ply absorb the hadronic contribution IIhad(g 2 ) into an 
alternate definition of the reduced form factors: 


F t (q 2 ) -A [1 - Hhadte 2 )]" 1 ^ 2 ) • (20) 


Several remarks are in order. First, we note that 


Eq. (161 is not a calculation of proton-vertex radiative 


corrections but rather a definition of Born form factors 
in the presence of radiative corrections. The definitions of 


the radii following from Eqs. (16) and (18) differ slightly 
from the definition of Maximon and Tjon [30) which in¬ 
cludes an additional contribution (there denoted dfj 1 ) in¬ 
volving a sticking-in-form-factors (SIFF) ansatz for the 
proton vertex. However, in many analyses (including, in 
particular, Ref. J), the additional contributions beyond 
those in Eq. (16) are anyways ignored. The convention 


(161 does not require the specification of a form factor 


model and is closely aligned with standard treatments of 
electron scattering. Let us further remark that this con¬ 
vention differs slightly from a convention commonly used 
in atomic physics applications [32j . However, the differ¬ 


ence, represented by the term a/{12Tvm 2 ) in Eq. (191, 


corresponds to a relative shift of ~ 2 x 10 -5 in te, well 
below current experimental sensitivities in either electron 
scattering or muonic hydrogen. 

Second, we remark that if hadronic vacuum polariza¬ 
tion is not removed explicitly before fitting then the re¬ 
sulting proton form factors should be interpreted with 
the alternate definition (20). With this definition, the 


fitted radius now corresponds to 


l r E,M\ * — r E,M + 6n; iad (o). 


( 21 ) 


A dispersive analysis of e + e —» hadrons data yields 


HLd(0) = -9.31(20) x 10- 3 GeV -2 . (22) 

This correction leads to a small shift, ~ 0.001 frn, in te, 
and a more careful error analysis does not appear to be 
warranted at the current level of precision. We note that 


' Alternatively, one could model the hadronic contributions by 
quark loop diagrams ED; however, strong interaction corrections 
are not controlled at small Q 2 < Aq CD . 
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Ref. [5] did not account for hadronic vacuum polariza¬ 
tion explicitly and hence implicitly employed the alter¬ 


nate definition (20). Experiments in the world data set 
used explicit models to account for hadronic vacuum po¬ 
larization or included uncertainties to account for the 
neglect of this correction. Hence, the extracted radii 
should differ slightly from the Mainz value according to 
the replacement (21), but the effect is well below the cur¬ 
rent experimental precision. 8 The treatment in Eq. (20) 
efficiently accounts for the effects of hadronic vacuum 
polarization on the radii in terms of the single number 
nfiad(O). When interpreting form factors at finite mo¬ 
mentum transfer, care must be taken to account for the 
q 2 dependence of hadronic vacuum polarization. 

Finally, we note that the analytic structures of the 
functions I\(p,p') in Eq. (16) and of ffhad(9 2 ) in Eq. (20) 
do not upset the assumptions going into the 2 expan¬ 
sion, since these functions are analytic outside of a cut 
at timelike q 2 > 4 m 2 for K(p,p') and at q 2 > 4m 2 for 

nhad(g 2 )- 


B. Two-photon exchange 

The two-photon exchange (TPE) contribution may be 
written 

- K(p, - k ) - K(p', -k') + K(p, k') 



TABLE I: Expansion coefficients for Eq. (251 in the SIFF 
TPE prescription of Ref. m- Note that 713 is determined by 

•Ei, 2 ( 0 ) = V.; "j/di- 



Fi 

f 2 

Til 

0.38676 

1.01650 

n 2 

0.53222 

-19.0246 

di 

3.29899 

0.40886 

d2 

0.45614 

2.94311 

d>3 

3.32682 

3.12550 


M 2 when applied to the Mainz data, which uses the same 
convention. For the world data, the Mo-Tsai convention 
is used, and so our calculation of _Ad^ [aTj yields a total 
M 2 contribution that differs by —0.4% to 0.1% at the 
cross section level compared to the consistent combina¬ 
tion. This small error is accounted for in the radiative 
correction uncertainties quoted for these measurements, 
and we have verified that such differences have an in¬ 
significant impact on the extracted radii. 

As a default, we employ the SIFF ansatz to estimate 
the TPE correction [35H35] . We have computed M 2 us¬ 
ing two form factor models. The first uses dipole FT, F 2 
form factors, 


Ei, 2 (<? 2 ) -> E 1]2 (0) 



(24) 


K(p',k) 


Mi+M 


MoTs 


E 




mi 


Tog 


'E+y/E 2 - ml 


E' 




log 


'e' + ^/e^i 


mt 


mi 


m e 


(23) 


x log -^2 Mi + Af“ aTj 


The A'(pi,p 2 ) functions are defined above in Eq. ( jll 
As indicated in Eq. (23), two conventions exist in the 


literature for isolating an IR finite TPE contribution: 
A4| loTs gg (Mo-Tsai) and Adf aTj (3D] (Maximon-Tjon). 
As long as the full correction M 2 is applied to the data, 
the results are independent of the convention used to 
separate IR divergent and IR finite contributions. Our 
hadronic model for the finite contribution is based on 
the Maximon-Tjon convention and so yields the complete 


with a value A 2 = 0.71 GeV 2 . The second model repre¬ 
sents Fi, F 2 as a sum of monopoles, 

< 25 > 

i=1 3 9 

To compare with previous results in the literature (35], 
we consider in particular the case N = 3 with parameter 
values rij , dj given in Table IT] We compare these models 
for M 2 to results with vanishing finite TPE correction in 
the Maximon-Tjon convention, 

A4^ aTj (no TPE) = 0, (26) 

and to results setting A4^ aTJ /AIi equal to the complete, 
“Feshbach” [35], result for M 2 /Mi in the rn p —► 00 
limit 9 : 

/ s j n l \ 

A4^ aTj (Feshbach) = I 1 +7ra- )Mi . (27) 

V 1 + sin |) 


An explicit correction for hadronic vacuum polarization is typ¬ 
ically applied in atomic physics analyses. This is the case in 
particular for the CREMA analysis of muonic hydrogen [3]- To 
avoid a double counting, the shift HD should therefore in prin¬ 
ciple be applied to electron scattering extractions that absorb 
this correction into the definition of the radius, before input or 
comparison to atomic physics extractions. 


9 An imaginary part in Eq. ( |27| is ignored since it affects the 
cross section only at relative order a 2 . For definiteness, we have 
expressed the result in terms of ( E , 0) instead of the variables 
(E, Q 2 ) before taking the m p —> oo limit, to match the expres¬ 
sion used in Ref. [9]. 
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C. Soft bremsstrahlung 


E. Summary of experimental implementations 


The soft bremsstrahlung contribution to the cross sec¬ 
tion is 


t^brem. — 


k' 


d 3 e 


V£ 2 + A 2 


\£\<(E/E')AE 


k'-t 


1 ’ (28) 


P 


V 


where A E is the accepted energy cut interval for the final 
state electron in the lab frame. This integral may be 
evaluated analytically [SD] ■ After cancellation of infrared 
divergences, the differential cross section including first- 
order real and virtual radiative effects may be written 


dcr = (da) o(l + (5). (29) 

Here, (da )o is the cross section ([5|, expressed in terms of 
Born form factors, and S is a finite correction depending 
on kinematic variables that accounts for vertex, vacuum 
polarization, and TPE radiative corrections. 


D. Large log resummation 


The preceding subsections (Secs. IV A TV C) summa¬ 
rize a complete treatment of first-order radiative correc¬ 
tions. The hadronic input, apart from the form factors 
to be determined, consists of a TPE model for M 2 in 
Eq. (231 and the number fth ad (0) (the latter impacts 
the radius at a level below current uncertainties). How¬ 
ever, we wish to describe scattering data with momen¬ 
tum transfers as large as Q 2 ~ lGeV 2 . In this regime, 
large logarithms from electron radiative corrections cause 
a poor convergence, or even breakdown, of the naive per¬ 
turbation theory, since 


«, 2 Q' 

-log — 
7 r m, 


0.5. 


Q 2 ~l GeV 2 


(30) 


Thus, first-order radiative corrections are insufficient for 
percent-level accuracy. 

When Q ~ E ~ E' and m e ~ A E, the leading series 
of logarithms a n log 2n (Q 2 /m 2 ) are resummed by making 
in (29) the replacement, 


1 + <5 —> exp (5). 


(31) 


Two-loop corrections without logarithmic enhancement 
are below the relevant experimental precision. For defi¬ 
niteness, in our analysis of the Mainz data, we employ the 
prescription used in Ref. )9], exponentiating all first-order 
corrections in (31) except the finite TPE contributions. 


We return to a discussion of deficiencies in this treatment 
in Sec. I VII ( 1 31 


Apart from TPE and, in some cases, large log resum¬ 
mation, radiative corrections have already been applied 
to all of the cross sections we include in our fit, as part of 
the original analysis of the experiments. We will exam¬ 
ine the impact of different TPE prescriptions, with final 
results based on the SIFF sum of monopoles TPE cor¬ 
rection as in Eq. (25) and Table [TJ Possible deficiencies 
in radiative corrections are treated at the same level as 
experimental systematic errors. 

Consider first the Mainz data set. The Al Col¬ 
laboration’s data analysis applied radiative corrections 
based on the prescription of Refs. m l40l . as detailed 
in Ref. [5]- This includes TPE corrections using the 
Feshbach prescription |39] and the large log resumma¬ 
tion given in Eq. (31) above (excluding the finite TPE 
contribution from the exponentiation). In the analysis 
of correlated systematic uncertainties, the cutoff on the 
bremsstrahlung tail was varied, yielding rms cross section 
variations well below 0.1%. These variations in the cross 
sections were used to determine the impact of radiative 
correction uncertainties on the radius. No uncertainty 
was included in the cross sections for the TPE contribu¬ 
tion. 


Consider now the world data set. The world data come 
from many different experiments, and the details of the 
radiative corrections vary. However, they are all based on 
the general formalism of Mo and Tsai 1411 312] , with 
improvements and modifications added in later works, 
e.g., Refs. [34] [43]. Our compilation of world data comes 
from Ref. [Hj, along with additional low Q 2 and more 
recent cross section Hul - FTS] and recoil polarization or po¬ 
larized target measurements [13] [49U56] , with earlier po¬ 
larization transfer results [57H59] replaced by the results 
of final, updated analyses USED- Further details of ra¬ 
diative corrections, in particular for earlier experiments, 
are presented in Ref. [63). Our compilation includes the 
corrections applied to earlier measurements discussed in 
that work; furthermore, we include additional vacuum 
polarization terms; exclude small angle data, 6 < 20°, 
from Ref. [43]; and include separate normalization fac¬ 
tors for data taken with different detectors or under very 
different conditions [MHflCl] . 

The original publications of the experiments compris¬ 
ing the world data set did not apply TPE corrections, and 
different prescriptions were used to approximate Eq. ©. 
For the most part, these experiments quoted normaliza¬ 
tion and uncorrelated uncertainties of 0.5-1% each to ac¬ 
count for uncertainties in the radiative corrections ap¬ 
plied, dominated by uncertainty associated with TPE 
corrections. 10 In this case, we will apply TPE correc- 


10 While this turned out to be smaller than the size of TPE cor¬ 
rections in recent calculations [33 [38], it appears to be a signifi¬ 
cant overestimate of the residual uncertainty at lower Q 2 values, 



















tions similar to those applied in the previous global anal¬ 
ysis [26j , for which the errors assigned in previous exper¬ 
iments were taken to be sufficient to account for uncer¬ 
tainties after applying a hadronic calculation of the TPE 
corrections. Note that one experiment [55] did not in¬ 
clude uncertainties associated with these corrections and 
so had much smaller total uncertainties than other exper¬ 
iments. Following Ref. [25], we thus include an additional 
systematic uncertainty to the data of Ref. [55]: we in¬ 
crease the normalization uncertainty by 1% (to a total of 
1.5%) and add 0.5% in quadrature to the point-to-point 
uncertainty. 

In Sec. m we will include constraints on the form 
factor ratio Ge/Gm from polarization measurements. 
In the kinematic range considered, with Q 2 < lGeV 2 , 
the TPE correction, estimated from a simple hadronic 
model [35], is small compared to experimental errors. * 11 
Following Ref. [25], this model-dependent correction is 
thus omitted from the fits. We will find that the po¬ 
larization data points do not have a strong influence on 
the radius fits, and thus do not pursue a more detailed 
treatment of radiative corrections to these data points. 


V. UPDATED FIT OF THE MAINZ DATA SET 

In this section, we extract the charge and magnetic 
radii from the Mainz data set, retaining the original treat¬ 
ment of statistical and systematic uncertainties and cor¬ 
rection factors from Ref. [9] but incorporating our knowl¬ 
edge of the structure of the form factors as presented in 
Sec. |III| We first reproduce the Mainz polynomial and 
inverse polynomial fits and then provide an updated ex¬ 
traction using the bounded z expansion. To highlight 
differences in the theoretical treatment of the form fac¬ 
tors, we fit to the full data set (1422 points) and apply 
the Feshbach correction as the only TPE correction, as 
was done in the primary radius extraction from Ref. [9]. 
We then discuss the impact of moving from polynomial 
fits to fits using the bounded z expansion and comment 
on other attempts to extract the radii from the Mainz 
data. 

Note that, because of the way the data and uncer¬ 
tainties are parametrized for the Mainz data, the uncer¬ 
tainties from such a fit represent only part of the total 
uncertainty. Meaningful error estimates require the ex¬ 
amination of correlated effects arising, e.g., from exper¬ 
imental systematic errors and radiative corrections. In 
this section, we focus on how the improved form factor 


based on the consistency between low-Q 2 estimates of the cor¬ 
rections m- 

11 The hadronic model for TPE corrections m predicts a correc¬ 
tion not larger than ~ 0.5% over the full e range for Q 2 < 1 GeV 2 . 
Furthermore, the 41 data points with Q 2 < 1 GeV 2 are concen¬ 
trated at large epsilon, e « 0.7 =t 0.2, where the TPE correction 
model predicts a correction < 0.2%. 


parametrization modifies the extracted radii and fit un¬ 
certainties. Section lYll will include an examination of the 
corrections applied to the Mainz data and the treatment 
of systematic uncertainties presented in their analysis. 

We take the cross section data as provided in the Sup¬ 
plemental Material of Ref. [9], which includes the Fes¬ 
hbach correction for TPE, and scaling of the statistical 
uncertainties to account for unidentified systematic er¬ 
rors, as discussed in Sec. |VI| We allow the normalization 
parameters to float freely in the fit, in accordance with 
Ref. [9]. In addition to examining the full Mainz data 
set, we also provide results obtained by restricting to 
momentum transfers below a given ax . We use the x 2 
function 

( / \ 2 

2 \ a i ~ a i, fit/^i,fit) / 0 r>\ 

x„ = E--■ < 32 > 

■;=l 2 


Here, N a is the number of cross section points for a 
specified kinematic cut ax , a is the measured cross 
section (after accounting for radiative corrections), 5a 
is the (point-to-point, uncorrelated) uncertainty, ergt is 
the cross section calculated using the chosen form factor 
model, and r/g t is a product of normalization parameters 
for a given run (i.e., for data taken at a given choice of an¬ 
gle and energy). There are 31 normalization parameters 
in the complete data set. 

In our default fits, we enforce bounds on the form fac¬ 
tor parameters by a % 2 penalty, 



(33) 


where a^gt and bi t g t are the fit values of the coefficients 
for G p e and G P M , respectively, and |a fe | max and |6 fe | max are 
(Gaussian) bounds on the coefficients. For the polyno¬ 
mial and unbounded z expansion fits, |a.fc| max and |&fc| max 
are taken to be very large, acting simply as numeri¬ 
cal regulators in the fits (they are taken large enough 
such that fit results represent the infinite bound limit). 
For the bounded z expansion, Eq. (331 enforces a Gaus¬ 
sian, vs sharp cutoff, statistical prior on the form factor 
parameter space, 12 typically taken to be |cifc/ao| max = 
|frfc/&o|max = 5. A more detailed discussion of the depen¬ 
dence of fit results on form factor priors is postponed to 
Sec. IVTTCl 


A. Polynomial and inverse polynomial fits 

The radius central values, minimum y , and reduced 
X 2 are displayed in Table [n] for fits with form factors 
represented as polynomials in q 2 of degree 10, or as in¬ 
verse polynomials of degree 7. These results are very 


12 For a related discussion, see Ref. m ■ See also Ref. EH- 
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TABLE II: Results for fits using polynomials of degree 10 and 
inverse polynomials of degree 7 for the full (N a = 1422) A1 
MAMI data set. The reduced y 2 is calculated taking Ndot = 
N a -2k 

max Anorm with N noim =31. 


Fit type 

r E [fm] 

r M [fm] 

x 2 

Xred 

poly 10 

0.886 

0.794 

1561.6 

1.14 

inv poly 7 

0.886 

0.768 

1569.1 

1.14 


close, but not identical, to the corresponding results in 
Table IV and Fig. 20 of Ref. ]5]- We have compared our 
results to the output from the example fitting code pro¬ 
vided as part of the Supplemental Material for Ref. 0, 
finding agreement with the results of this code. For ex¬ 
ample, in the case of the polynomial of degree 10, the 
results of the example fitting code agree with our results 
in Table [TT} both having a minimum y 2 °f 1561.6, lower 
than the value 1563 quoted in Table IV of Ref. (Sj. 13 


B. Bounded z expansion fits 


Let us proceed to consider the implications of the 
bounded z expansion. Here we retain the identical data 
set as employed in Table [TT] For the default fit, we 
take to = 0, k max = 12, and a Gaussian bound of 
(, k |max = 5/. |max//q> = 5. The value fc max = 12 is 
large enough that the result does not change if k max is 
increased further. 


The results for this fit are displayed in Fig. [I] as a 
function of Q^ax- The extracted radii a nd y 2 values are 


provided for three Q^ ax values in Table III The quoted 


uncertainty includes only the statistical-type uncertain¬ 
ties, i.e., counting statistics and uncorrelated systematic 
uncertainties that are represented by the rescaling of the 
statistical errors in the Al data set. The uncertainty is 
obtained by varying the radius around the best-fit value, 
refitting the data while allowing all data set normaliza¬ 
tions to float, to map out the y 2 contour as a function of 
radius. The contours are typically symmetric and very 
nearly parabolic, and in the tables, we quote the aver¬ 
age of the change in radius that yields Ay 2 = 1 on 
the high and low sides of the central value. Note that 
the primary Al analysis of the Mainz data, identical 
except for the choice of the fitting function, yielded [5] 
r E = 0.879(5) sta t fm and r M = 0.777(13) stat fm, includ¬ 
ing only statistical uncertainties for comparison with our 
bounded 2 expansion results in Table |III[ 


13 More precisely, the fitting code returned a y 2 of 1561.60 and 
r M = 0.797fm. Evaluating our y 2 function with the correspond¬ 
ing parameters yielded an identical 1561.60. Using the same 
initialization conditions as the example fitting code, our mini¬ 
mization code independently returned a minimum y 2 of 1561.58 
and r ..y/ = 0.794 fm, as displayed in Table [TT] 



GLx [GeV 2 ] 

FIG. 1: Extracted electric (top panel) and magnetic (bot¬ 
tom panel) radii as functions of the kinematic cut ax on 
momentum transfer for the 1422 point Al MAMI data set, 
using the z expansion with to = 0, Gaussian priors with 
IUfc|max = |bfc|max/Mp = 5, fc max = 12. One-cr error bands 
are statistical only. 


TABLE III: Results from the fits in Fig. [I] for three values of 
Gmax- No- is the number of cross section points with Q 2 below 
Omax 5 and N n omi is the number of normalization parameters 
appearing in the data subset. 


Graax (GeV 2 ) 

r E (fm) 

r M (fm) 

Xmin 

N a 

N norm 

0.05 

0.873(18) 

1.071(114) 

479.4 

483 

13 

0.5 

0.905(10) 

0.749(28) 

1404.7 

1285 

29 

1 

0.920(9) 

0.743(25) 

1605.5 

1422 

31 


C. Discussion 

Let us remark on three aspects of the fits summarized 
in Table |TTT| First, we remark that the bounded z expan¬ 
sion fit to the entire 1422 point data set (Qm ax = 1 GeV 2 ) 
yields an electric radius significantly larger than the 
Mainz Al extraction [5]. Having analyzed identical data 
sets, this difference arises solely from requiring the form 
factors to lie within the class allowed by the bounded z 
expansion. The difference, 0.041 fm, is large compared to 
the Mainz estimated systematic uncertainty. The mag¬ 
netic radii exhibit a smaller difference, with our result 
0.034 fm below the Mainz extracted value. 

Second, the extracted radii have significant depen¬ 
dence on Q 2 nax . For example, r E = 0.873(18) sta t fm 
with Qmax = 0-05 GeV 2 vs r E = 0.920(9) sta t fm with 
Qmax = 1 GeV 2 . The difference, 0.047 fm, is again 
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large compared to the quoted uncertainties. Further¬ 
more, there is a non-negligible variation of the te central 
value as <5^ ax is increased above 0.5 GeV 2 , even though 
the region below 0.5 GeV 2 includes more than 90% of the 
data points, and (as illustrated below in Fig. 10) the data 
above 0.5 GeV 2 do not significantly impact the radius un¬ 
certainty. In fits with unbounded parameters, it is not 
surprising that the extracted radius is sensitive to higher- 
Q 2 data because the radius may change to provide a bet¬ 
ter fit to fluctuations in the data that are accommodated 
by arbitrarily large parameter values. This behavior is 
unexpected in fits with bounded parameters. Thus, it is 
surprising that the small amount of higher-Q 2 data has 
such a significant impact on the extracted radius. The 
dependence on Q^ ax suggests a possible tension between 
the lower- and higher-Q 2 data. 


Third, taking at face value the complete 1422 point 
data set and error assignments, the resulting electric ra¬ 
dius is te = 0.920(9) sta t(6)other fm, where for the mo¬ 
ment we simply take the A1 evaluation of other contri¬ 
butions to the uncertainty. 14 This result is 7a above the 
muonic hydrogen value, te = 0.84087(30) fm |J. It is 
also in tension with the results extracted from hydrogen 
spectroscopy, te = 0.8758(77) fm P, and with a previ¬ 
ous global analysis [13] of world electron-proton scatter¬ 
ing data which yielded r E = 0.875(10). The magnetic 
radius value of tm = 0.743(25) s t a t(10) o ther is almost 4 a 
from the value = 0.867(20) fm from the global anal¬ 
ysis in Ref. [13] . However, we note that recent global 
analyses Ham, use different representations of the form 
factors compared to the bounded z expansion used in 
Table III In Sec. VIIB below, we will perform our own 
analysis of the world data for a consistent comparison 
with the analysis of the Mainz data. 


1750 

1700 

1650 

1600 



FIG. 2: Total y 2 (top panel) and extracted electric (middle 
panel) and magnetic (bottom panel) radii as functions of fc max 
for the 1422 point A1 MAMI data set, using the z expansion 
with to = 0, Gaussian priors with |a fc | max = |6fc| max //r p = 5. 
One-<T error bands are statistical only. 


D. Further tests related to the z expansion 


1. Dependence on fc max 


Simply replacing the fit functions employed in Ref. [5] 
with the z expansion does not resolve the discrepancy 
with muonic hydrogen results. In fact, the result is a 
larger difference with muonic hydrogen, as well as a ten¬ 
sion with previous extractions from world electron-proton 
scattering data. In addition, the results show an unex¬ 
pected dependence on the Q 2 range of data included in 
the fit. In the following Sec. |VI[ we consider in detail a 
range of sources of the systematic error before presenting 
best values for the radii. 


14 The error (6) ot her results from the quadrature sum for the errors 
(4)syst (2) model (4) group presented in Ref. [9]. These errors were 
added in quadrature in Ref. [9], but it has more recently been 
advocated m to add the final error linearly to the quadrature 
sum of the first two, resulting in ~ (8.5) 0 ther- 


In the bounded z expansion, we may estimate the max¬ 
imum power of z which can impact the data at a given 
level when the expansion coefficients a n are order unity. 
Setting the upper limit of the contribution at the level of 
~ 0.5% implies fc max ~ 10 should be sufficient for t 0 = 0, 
Qma.x = 1 GeV 2 . Figure [ 2 ] shows the y 2 values and radii 
extracted as a function of fc max for the bounded z ex¬ 
pansion fit to the full Mainz data set. The rightmost 
points at fc max = 12 correspond to the rightmost points 
in Fig. [T] and to the final row of Table [TTT] In accordance 
with our power counting estimate, the minimum y 2 value 
and extracted radii have stabilized by & max = 10. For 
definiteness, we choose fc max = 12 for all of our bounded 
z expansion fits. While this is significantly more param¬ 
eters than required for fits to smaller values of Q 2 lax , the 
bounds on the fit parameters prevent the problem of ra¬ 
dius instability due to overfitting of noise in the higher-Q 2 
data [IT] . 
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2. Unbounded z expansion fits 




FIG. 3: x 2 (top panel, solid black line) and extracted elec¬ 
tric (middle panel, solid black lines) and magnetic (bottom 
panel, solid black lines) radii with la statistical error bands 
as functions of fc max for the unbounded z expansion fit with 
to = 0 to the 1422 point A1 MAMI data set (with floating 
normalization). The dashed blue lines show the x' 2 and 1 a 
error bands from the bounded fit in Fig. [2] for comparison. 


The bounded z expansion (the formal fc max —> oo 
limit with bounded coefficients) is a particularly well- 
motivated implementation of form factor priors. A dif¬ 
ferent but common choice of priors corresponds to set¬ 
ting dfc = 0 for all coefficients beyond a given order 
k > fc max , with the remaining coefficients unconstrained, 
and —oo < cifc < oo for k < fc max . We perform some 
illustrative fits with this modified choice of priors in or¬ 
der to separate the impact of applying bounds from the 
impact of changing from polynomial or inverse polyno¬ 
mial functions to the z expansion. We again fit the 
1422 point Al data set using the same rescaled errors 
and Feshbach TPE correction as in Fig. [Tj but now set 
| tlk | max 5 |^fc|max t OO in Eq. (33). 

In the limit of large fc max , the true form factors are 
guaranteed to lie in the space of curves described by the 
unbounded z expansion. However, many badly behaved 
form factors (in particular, form factors in conflict with 
predictions of QCD, as discussed in Sec. Ill) also lie in 


this space of curves, and fits without constraints on the 
coefficients lose predictive power at large fc max . 

Figure [3] shows results for unbounded fits with float¬ 
ing normalizations. The minimum y 2 value continues to 
decrease significantly as parameters are added through 
fc m ax = 10. Quantitatively reliable radius estimates are 
difficult to obtain from such a fit; for small fc max , omit¬ 
ted terms in the form factor expansion can introduce a 
potentially large, but difficult to quantify, bias in the fit¬ 
ted radii [71], while for large fc maX) the uncertainties grow 
rapidly. 


3. Fixed-normalization fits 




FIG. 4: x 2 (top panel) and extracted electric (middle panel) 
and magnetic (bottom panel) radii with la statistical error 
bands as functions of fc max for the unbounded a expansion fit 
with to = 0 to the 1422 point Al MAMI data set with fixed 
normalization parameters. 


As noted earlier, the manner in which uncorrelated sys¬ 
tematic uncertainties are treated in the Mainz data set 
is only complete when the normalization parameters are 
allowed to vary and when correlated systematic uncer¬ 
tainties are estimated separately [9]. Thus, fits which fix 
the normalization of the data based on the default Mainz 
fit will underestimate uncertainties and potentially yield 
different values for the radii if the fit is performed with a 
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different functional form than that used to determine the 
normalization parameters. Nonetheless, such fits have 
been performed 0(72], and so we provide a comparison 
of unbounded fits with and without floating normaliza¬ 
tion factors. 

The result for fixed normalization factors is displayed 
in Fig. [4j Comparing to fits with floating normaliza¬ 
tions, one can see that the uncertainties are significantly 
smaller in the case of fixed normalizations, with the fit 
uncertainties on rg underestimated by a factor > 5. 15 

Even ignoring the artificially small uncertainties that 
arise when neglecting the normalization uncertainty of 
the data sets, it is not clear that there is any value of /c max 
for which the fit provides a sufficiently precise description 
of the data while still providing meaningful uncertainties 
on the charge radius. For the magnetic radius, the results 
are even less clear, with only an upper limit on the radius 
possible for fc max > 9. We note that for large fc max the 
central values for the fits displayed in Figs. [3] and [4] 
require very large coefficient values, in violent conflict 
with order unity predictions of QCD. lh 


VI. SYSTEMATIC STUDIES FOR THE MAINZ 
DATA SET 


Taking the data and error prescriptions of Ref. [9[ at 
face value, we have found radius extractions in tension 
with each other for fits of different functions to the same 
data set (Table [II] compared to the last line of Table III) 
and for fits of the same function to subsets of the same 
data set (Fig. [I] and Table |HI[ ). 

The first observation indicates the strong dependence 
of the extracted radius on the specification of physical 
form factors, as implemented by the bounded z expan¬ 
sion. We focus solely on this class of form factors in the 
following. Since the systematic error analysis of Ref. [5] 
relied on rescaling statistical errors based on fits to a par¬ 
ticular class of form factor models, we also revisit this 
analysis using the bounded z expansion. 

The second observation, that fits to data subsets are 
in tension with fits to the entire data set, indicates the 
possibility of an underestimated systematic error. We 
investigate below a range of correlated errors and their 
potential impact on radius extractions. 

The analysis of the Mainz data by the Al Collabora¬ 
tion decomposed the uncertainties into several different 
contributions. Let us briefly review this decomposition. 
The only uncertainty applied directly to the quoted cross 
sections is called the statistical uncertainty. This is a 


15 In detail, for fc m ax = 9(10), the radius uncertainty is Ste = 
0.011(0.014) for fixed normalizations, compared to 8te = 
0.053(0.096) for floating normalizations. 

16 For example, at fc max = 9, requiring the central value for 77 ,; in 
Fig. a to lie within the 1 rr envelope of a bounded z expansion 
requires coefficient bound > 10 4 . 


combination of counting statistics and systematic uncer¬ 
tainties which are taken to be uncorrelated between dif¬ 
ferent data points and normally distributed and is thus 
treated in the same way as counting statistics. We refer 
to the systematic uncertainties that are independent for 
different data points as the uncorrelated systematic un¬ 
certainties. The Al Collaboration includes these uncorre¬ 
lated systematic uncertainties by introducing a rescaling 
factor on the counting statistics, with a procedure to ex¬ 
tract these scaling factors which we summarize below in 
Sec. I VIB II 

Normalization uncertainties for each data subset in the 
experiment are accounted for in the extraction of the ra¬ 
dius by allowing the 31 normalization factors correspond¬ 
ing to different configurations to float freely when fitting 
the form factors. The Al analysis of Ref. [S] suggests 
an uncertainty of 3.5-5% on the normalization factors, 
but no constraints were included in the fits. Because the 
cross sections are quoted after the determination of the 
normalization factors in their fit, any information on the 
initial normalizations is lost, and it is no longer possi¬ 
ble for us to make use of even the limited precision with 
which these normalization factors were constrained. 

Finally, the Al Collaboration gives a procedure for es¬ 
timating the impact of correlated systematic uncertain¬ 
ties on their data. These are corrections which poten¬ 
tially have a strong kinematic dependence, and their im¬ 
pact will not necessarily decrease when one includes a 
large number of measurements. Such corrections must 
be treated independently from the statistical and uncor¬ 
related systematic errors. The Mainz treatment of these 
uncertainties and our examination of other possible cor¬ 
related effects are included in Sec. I VI Cl 

In the remainder of this section, we introduce the 
following three modifications to the analysis. First, in 
Sec. |VI A[ we study the impact of different TPE cor¬ 
rection models on radius extractions, choosing the SIFF 
sum of monopoles ansatz as the default in the remain¬ 
ing fits. Second, after identifying in Sec. VIB 1 poten¬ 
tial shortcomings in the rescaling of statistical errors, in 
Sec. |VI B 2[ we rebin data taken at identical kinematic 
settings in order to incorporate in Sec. VIB 3 uncorre¬ 


lated systematic errors which do not scale with statistics. 
Lastly, in Sec. | VIC] we consider a range of correlated sys¬ 
tematic errors consistent with the experimental precision 
achieved in Ref. |J)]- 


A. TPE model dependence 


The variation of the extracted radii with different TPE 
models is illustrated in Fig. [5] where the extracted te and 
r M central values are plotted vs ax - The black curve 
is identical to the central curve in Fig. [I] The remaining 
results are obtained by repeating the fit of Sec.[V]after re¬ 
moving the Feshbach TPE correction, Eq. ( [27] ) , and then 
applying the SIFF TPE result [using dipole form factors, 
Eq. (24), or those from Blunden et al., Eq. (25)] or ap- 
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Qls* [GeV 2 ] 


FIG. 5: Extracted electric (top panel) and magnetic (bottom 
panel) radii as functions of the kinematic cut QJkax on momen¬ 
tum transfer for several TPE models, as discussed in the text: 
no correction (red, dotted), Feshbach correction (black, solid), 
SIFF dipole (green, dot-dashed), and SIFF sum of monopoles 
(blue, dashed). There is a negligible difference between the 
SIFF choices of the dipole and the sum of monopoles. Fits are 
to the 1422 point A1 MAMI data set, using the 3 expansion 
with to = 0, Gaussian priors with |a fc | max = |&fc[ ma x//+> = 5, 

max — 12. 


mal limit of infinite proton mass and is independent of 
the proton structure. The exact result for arbitrary kine¬ 
matics for a pointlike proton [35] yields a correction that 
grows with Q 2 , approximately doubling the correction 
between Q 2 = 0 and 1 GeV 2 . However, calculations us¬ 
ing either hadronic }35j or partonic m models to account 
for proton structure indicate that the correction does not 
grow with increasing Q 2 but instead becomes smaller and 
then changes sign. This is the behavior required to ex¬ 
plain the difference between the Rosenbluth and polar¬ 
ization measurements of [i v Ge/Gm for the proton [ 23 ] 
and has been recently confirmed for Q 2 « 1 T.5 GeV 2 by 
comparisons of positron and electron scattering from the 
proton [731 EH- 

There is a significant difference in the charge radius be¬ 
tween the case of no TPE corrections and either the Fes¬ 
hbach or SIFF corrections. However, there is a relatively 
small difference between Feshbach and SIFF, suggesting 
that the infinite proton mass limit provides a significant 
part of the correction for r^. For the magnetic radius, 
there is a large difference between all three approaches. 
For both the charge and magnetic radii, there is little 
sensitivity to the choice of form factors included in the 
SIFF calculation. We collect in Table m the deviations 
of the extracted radius using different models in place of 
the Feshbach correction. In all subsequent fits we em¬ 
ploy the SIFF ansatz, using for definiteness the sum of 
monopoles in Table [T| as our default TPE model. The 
uncertainty associated with TPE corrections will be in¬ 
corporated into the evaluation of correlated systematic 
uncertainties in Sec. I VI Cl 


TABLE IV: Change in the extracted charge and magnetic 
radii for three different TPE corrections, relative to the Fes¬ 
hbach correction applied in the Mainz analysis. Results are 
for the fit with ax = 0.05, 0.5, 1 GeV 2 in Fig. [H] 


Qmax (GeV 2 ) 

Model 

A r E (fm) 

A r M (fm) 

0.05 

Feshbach 

= 0 

= 0 


SIFF dipole 

-0.004 

+0.022 


SIFF Blunden 

-0.004 

+0.025 


No TPE 

-0.023 

-0.028 

0.5 

Feshbach 

= 0 

= 0 


SIFF dipole 

-0.003 

+0.036 


SIFF Blunden 

-0.002 

+0.034 


No TPE 

-0.017 

-0.026 

1 

Feshbach 

= 0 

= 0 


SIFF dipole 

-0.003 

+0.038 


SIFF Blunden 

-0.002 

+0.037 


No TPE 

-0.016 

-0.026 


plying no finite TPE correction [in the Maximon-Tjon 
convention, Eq. (26)]. As the plot illustrates, expressed 
as a difference relative to the Feshbach correction, the re¬ 
sults have mild Qmax dependence. Numerical values for 


Quia* = 0.05, 0.5, 1 GeV 2 are given in Table IV 


The Feshbach correction is the exact result in the for- 


B. Uncorrelated systematic uncertainties 

1. Summary of the Mainz A1 approach 

To estimate the uncorrelated systematic uncertainties, 
the A1 Collaboration performed a fit to the entire 1422 
point data set using a default form factor model (an 
eight-parameter cubic spline model for each of Ge and 
Gm)- The data were then grouped according to the beam 
energy and the spectrometer used in the measurement. 
For each data group, the uncorrelated systematic uncer¬ 
tainties were taken from examination of the distribution 
of the differences between measured and fit cross sections, 
scaled by the uncertainty from counting statistics. (If the 
counting statistics fully represented the uncorrelated un¬ 
certainties, then this should be a Gaussian distribution 
with width one.) This distribution was fit with a Gaus¬ 
sian, the width of which was then taken as the scaling 
factor applied to the counting statistics to determine the 
combined statistical and systematic uncorrelated uncer¬ 
tainties. The scaling factors obtained in this way vary 
from 1.070 to 2.283, as given in the Supplemental Mate¬ 
rial of Ref. 15]. 

This rescaling procedure is meant to yield a reduced 
X 2 close to unity when the data are compared to the 
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original fit. However, because the Gaussian fit may un¬ 
derestimate the impact of outliers and the scaling of the 
uncertainties changes the relative weighting of the differ¬ 
ent data sets, the fit to the data set with updated uncer¬ 
tainties yields a reduced y 2 somewhat larger than unity: 
X'red ~ 1-15 for the entire data set. This suggests that 
the quoted systematics are somewhat underestimated. 

The rescaled statistical errors represent the minimum 
additional uncertainty necessary to account for random 
scatter around a global fit to the data. Any correlated 
effects will not be included in the extracted uncertainties. 
For example, the Al analysis of Ref. [9] does not include 
any uncertainty associated with the error in the measure¬ 
ment of the beam energy or offsets in the spectrometer 
angles. Such kinematic offsets would yield correlated er¬ 
rors in the cross sections which would not be captured 
by this procedure. 

The Al rescaling procedure yields systematic uncer¬ 
tainties which depend on the form factor model used 
in the fit. We performed a similar analysis using our 
bounded z expansion and using the y 2 value for each 
data subset relative to the fit as the square of the scal¬ 
ing factor. This procedure yielded similar scaling factors, 
larger by 6% on average compared to the Al procedure 
(thus yielding a value of y 2 ed closer to unity), with a 
typical scatter around the small average offset of approx¬ 
imately 10%. Thus, the form factor model dependence of 
this rescaling procedure is small, though not negligible, 
and related more to the difference in our procedure than 
to the change in the fitting function. 

The rescaling procedure accounts for undetermined 
systematic errors that are assumed to be uncorrelated 
and to scale with the statistical counting errors. How¬ 
ever, it is not clear that all the uncorrelated systematics 
should scale with the statistical uncertainties. If one as¬ 
sumes that the statistical and uncorrelated systematic 
uncertainties add in quadrature, then one can compare 
the original (unsealed) and rescaled statistical uncertain¬ 
ties to extract the effective uncorrelated systematic un¬ 
certainty used in the Al procedure. This inferred un¬ 
certainty can be extremely small, as low as 0.05%, but 
varies with kinematics and with the counting statistics, 
attaining values up to 2%. Because the full set of 1422 
data points includes many instances of repeated mea¬ 
surements at identical kinematics, this procedure implies 
even greater reduction in the systematic uncertainty asso¬ 
ciated with each independent kinematic point, with val¬ 
ues as low as 0.02%, even though it is experimentally 
difficult to constrain uncertainties at that level. 

To address these concerns, we present a modified treat¬ 
ment of the uncorrelated systematic uncertainties where 
a fixed uncorrelated systematic error is added to all 
points to account for unknown drifts or corrections. 


2. Rebinning studies 

As noted above, the 1422 data points include many re¬ 
peated measurements at the same conditions. One would 
expect that many potential systematic errors would be 
identical for these points, e.g. time-dependent efficien¬ 
cies, rate-dependent corrections, or uncertainties in the 
beam energy or spectrometer angle settings. Adding a 
fixed systematic to every one of the 1422 data points 
would underestimate the systematic uncertainty for data 
points with many repeated measurements. Therefore, we 
begin by combining data points taken with identical con¬ 
ditions, reducing the data set from 1422 data points to 
658 independent cross section measurements. 

We group (i.e., rebin) all data taken at identical kine¬ 
matic settings, using only the uncertainties from counting 
statistics (i.e. removing the Al scaling factor). We tested 
our assumption that the points within the groups of re¬ 
peated measurements were consistent within statistics by 
looking at the y 2 values and confidence levels for every 
set of the rebinned data. There were 407 settings with 
multiple runs taken under identical conditions, and the 
confidence-level distribution for these sets is consistent 
with a uniform distribution between 0 and 100% except 
for a handful of outliers below 1% confidence level, indi¬ 
cating a nonstatistical scatter of the points being com¬ 
bined. Most of these outliers involved scatter at the 0.1 
0 .2% level, which is contained within the systematic un¬ 
certainty we will add to achieve y 2 ed near unity. One 
setting—Abeam = 315 MeV, 9 = 30.01°, spectrometer 
C—had a single measurement that deviated by ~1.5% 
from the two other measurements at that setting, while 
the statistical uncertainties were approximately 0.15%. 
We excluded this set of points, yielding a total of 657 in¬ 
dependent cross section measurements when fitting the 
rebinned data. 

We remark that normalization parameter 14 appears 
for only one point in the rebinned 657 point data set (and 
two points in the original 1422 point data set). Since 
the normalization parameters have to be allowed to float 
freely in the fit, this data point has no impact, but for 
definiteness, it is retained. Note that with the Mainz pro¬ 
cedure of applying scaling factors to the counting statis¬ 
tics the rebinning has no impact on the fit, and the only 
change at this point would be due to the exclusion of 
the one point after rebinning. However, when applying a 
more conventional constant uncorrelated systematic un¬ 
certainty to all points, the uncertainty is best applied to 
the rebinned data points. 


3. Uncorrelated systematics for rebinned data 

With the rebinned data set in hand, we proceed to 
investigate the inclusion of an uncorrelated systematic 
error that does not scale with statistics. We add a fixed 
systematic uncertainty to every data point and perform 
the bounded z expansion fit with Qmax = 1 GeV 2 and 
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TABLE V: Number of data points, reduced y 2 , an< i confi¬ 
dence level for each combination of spectrometer (A, B, or 
C) and beam energy (in MeV) of the rebinned A1 MAMI 
data set. Columns 4 and 5 give the results after the inclusion 
of a uniform 0.25% uncorrelated systematic; columns 6 and 
7 give the results after the inclusion of the final 0.3%-0.4% 
uncorrelated systematic. See the text for details. 


Spec. 

Beam 

N a 

Xred 

CL (%) 

Xred 

CL (%) 

A 

180 

29 

0.59 

96.1 

0.46 

99.4 


315 

23 

0.54 

96.4 

0.44 

99.1 


450 

25 

1.52 

4.8 

1.00 

46.7 


585 

28 

1.54 

3.4 

1.03 

42.8 


720 

29 

1.05 

39.9 

0.87 

66.4 


855 

21 

0.92 

56.8 

0.77 

76.0 

B 

180 

61 

0.85 

79.8 

0.65 

98.3 


315 

46 

1.05 

38.5 

0.76 

88.5 


450 

68 

0.90 

71.7 

0.67 

98.2 


585 

60 

0.61 

99.2 

0.50 

99.96 


720 

57 

1.29 

6.9 

0.97 

53.7 


855 

66 

1.88 

0.002 

1.15 

19.6 

C 

180 

24 

0.88 

63.3 

0.68 

88.0 


315 

24 

1.16 

27.2 

0.78 

76.8 


450 

25 

1.53 

4.3 

1.08 

35.9 


585 

18 

0.83 

66.3 

0.65 

86.4 


720 

32 

1.11 

30.2 

0.90 

62.3 


855 

21 

0.79 

73.7 

0.62 

90.5 


our default form factor scheme, to = 0, k max = 12, and 
Gaussian bound |afc|max = |6fc|max/ = 5. We varied the 
systematic uncertainty until we found a reduced y 2 value 
close to unity. This required a systematic uncertainty of 
approximately 0.3%. We then examined the y 2 contribu¬ 
tion from each of the 18 energy-spectrometer combina¬ 
tions to see if any of them had anomalously large or small 
Xied va l ues - While the spread of y 2 ed values was signifi¬ 
cant, many data subsets have a relatively small number of 
points, and the only subset which was an extreme outlier 
was the data from spectrometer B at E beam = 855 MeV. 
We chose to increase the systematic uncertainty on this 
data subset to 0.4%, while keeping 0.3% for all other data 
subsets. The reduced y 2 and confidence levels for each 
data sub set are displayed in Table [v] The total y 2 is 
520.4 for 657 points, which might suggest that 0.3% is 
a slight overestimate of the uncorrelated systematic, but 
it is a small effect, with a 0.25% correction yielding a 
reduced y 2 above 1 by a similar amount. 

Table IVII shows the radius fit results for the rebinned 
Mainz data with the statistical scaling factors from the 
original analysis replaced by the constant 0.3% system¬ 
atic uncertainty (0.4% for spectrometer B at 855 MeV 
beam energy). 

This procedure introduces enough uncertainty to ac¬ 
count for random scatter of the points around the best-fit 
curve. However, any errors that are correlated between 
multiple points will bias the fit, and will not be fully 
reflected in this procedure, making the resulting uncer¬ 
tainty estimate more of a lower limit. While the impact 


TABLE VI: Results for fitting of the 657 point rebinned A1 
MAMI data set with 0.3%-0.4% uncorrelated systematic un¬ 
certainties at three values of ax using the z expansion with 
to = 0, Gaussian priors with |afc| ma x = |&fc|max/7tp = 5, 
fcmax = 12. N a is the number of cross section points with 
Q 2 < Qmax, and N n orm is the number of normalization pa¬ 
rameters appearing in the data subset. 


QLx (GeV 2 ) 

te (fm) 

r M (fm) 

Xmin 

N a N nolm 

0.05 

0.856(27) 

1.11(14) 

110.5 

176 13 

0.5 

0.895(14) 

0.777(34) 

442.0 

568 29 

1 

0.908(13) 

0.767(33) 

520.4 

657 31 


of correlated uncertainties will be examined separately, 
these rely on specific models for kinematic dependences 
of any additional errors. The inclusion of an even larger 
uncorrelated uncertainty would allow the data to account 
for a range of correlated errors, but the reduced y 2 would 
end up significantly smaller than unity. For illustration, 
Table VII shows the results where we apply a 0.5% un¬ 
correlated systematic uncertainty to every data point, 
instead of the 0.3%-0.4% uncertainties in the previous 
fit. 


TABLE VII: Same as Table |VT| but with 0.5% uncorrelated 
systematic uncertainty. 


Qmax (GeV 2 ) 

rs (fm) 

r M (fm) 

Xmin 

N a 

N norm 

0.05 

0.861(35) 

1.05(18) 

48.7 

176 

13 

0.5 

0.891(18) 

0.768(43) 

211.5 

568 

29 

1 

0.901(17) 

0.758(42) 

250.3 

657 

31 


C. Correlated systematic uncertainties 

We now consider systematic errors that do not scale 
with statistical errors but that are also correlated across 
data points. We begin by examining the procedure of 
Ref. j5j. We then examine modified evaluations of the 
radius uncertainty associated with the correlated system¬ 
atic uncertainties. 


1. Summary of the Mainz A1 approach 

In the A1 MAMI data set, each cross section is accom¬ 
panied by two factors to account for systematic uncer¬ 
tainties. The first is due to the bremsstrahlung energy 
cut and is estimated by varying the cut. The second is 
meant to account for efficiency changes, normalization 
drifts, variations in spectrometer acceptance, and back¬ 
ground misestimations. This second class of systematics 
is investigated by applying a kinematic-dependent cor¬ 
rection to the data. The complete data set is refit after 
multiplying or dividing the individual cross section ratios 
by the corresponding factor for either the energy cut or 
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TABLE VIII: Results for changes in the radii under increases 
(upper value for each Qmax) or decreases (lower value) in the 
energy loss cut. Fits are for the 657 point rebinned A1 MAMI 
data set with 0.3%-0.4% uncorrelated systematic uncertain¬ 
ties at three values of ax using the z expansion with to = 0, 
Gaussian priors with |a fc | max = |&fc|max/M P = 5, fc max = 12. 


Ql ax (GeV 1 2 ) 

A r E (fm) 

A r M (fm) 

0.05 

-0.001 

+0.023 


-0.005 

0.000 

0.5 

+0.003 

+0.003 


-0.003 

+0.003 

1 

+0.003 

+0.009 


-0.002 

0.000 


correlated systematic error, and the largest difference in 
radius obtained from multiplying or dividing is taken as 
the uncertainty. The total systematic uncertainty is then 
obtained by summing in quadrature: 

Ar syst = \J (Ar Ecu t) 2 + (Ar corr ) 2 . (34) 


The stated cross section uncertainties associated with 
the variation in the energy cut are small, with a rms 
variation of 0.08%. These mainly introduce an additional 
scatter into the cross sections but have little impact on 
the radius central values. For the entire data set, this 
translates to an uncertainty in te of 0.003 frn and in tm 


of 0.009 fm. Explicit results are given in Table VIII 


In the Al analysis, the kinematic-dependent correlated 
systematic is assumed to depend linearly on the scat¬ 
tering angle [cf. Eq. (35) below, with x = 0], with a 
variation of approximately 0.2% between the minimum 
and maximum angles for each energy-spectrometer com¬ 
bination, except the 855 MeV data with spectrometer 
C (covering large angles), for which the variation is ap¬ 
proximately 0.5%. 17 We perform a more comprehensive 
study of correlated systematics below. 


2. Sensitivity to size or kinematic dependence 

The correlated systematics mentioned above could rep¬ 
resent either experimental or theoretical uncertainties. 
For example, they could be associated with radiative cor¬ 
rections (beyond the energy cutoff variation), background 
subtraction Mi potential offsets in the absolute beam 
energy or angle calibration, etc. The impact of such un¬ 
certainties on the cross sections is difficult to constrain 
below the 0.5% level, but because of the floating normal¬ 
izations of the different datas ets, these correlated sys¬ 
tematic uncertainties only need to account for the varia¬ 
tion within a specific normalization subset. 


1 ' These values are deduced from the appropriate column of the 

tabulated data set in the Supplemental Material of Ref. [9]. 


While some sources of correlated corrections may be 
well approximated by a correction that is linear in the 
scattering angle over a single energy-spectrometer set¬ 
ting, this is not the only possible kinematic dependence, 
and effects may be relevant over larger or smaller subsets 
of data, or may be more important for one spectrome¬ 
ter. Thus, we examine the impact of different prescrip¬ 
tions for applying the correlated systematics. We take 
a 0.5% variation in the systematic correction, but vary 
the functional form used to go from the minimum to the 
maximum kinematic settings within data subsets, and 
we vary how the full experiment is broken down. For the 
latter, we examine three cases: 0.5% variation over the 
range of angles for each spectrometer-energy combina¬ 
tion (as done in the Al analysis, with 18 separate angu¬ 
lar ranges), 0.5% variation over the full kinematic range 
for each spectrometer (with three separate ranges), and 
0.5% variation for each of the 34 normalization subsets. 

We examine eight different approaches to varying the 
kinematic dependence of the systematic correction over 
a given data subset. We multiply and divide the cross 
sections and uncertainties by the factor 


1 T dcorr — 1 




(35) 


where a = 0.005 and x is a kinematic variable. We take 
the variable x to be proportional or inversely proportional 
to 9, Q 2 , or E' , or to be proportional to e or 1/ sin 4 (6 ) /2). 
Note that for a given energy the correction goes from 
zero at one extreme of the angular range for the data 
subset to 0.5% at the other extreme; these different cor¬ 
rections only modify the interpolation to intermediate 
angles. These illustrative functional forms can be moti¬ 
vated from specific sources, including kinematic offsets, 
rate-dependent effects, or simplified models of radiative 
corrections. However, the exact magnitude and precise 
functional form cannot be fully determined without fur¬ 
ther input. 

Taking the correction to be linear in the scattering 
angle, x = 9, and applied to each of the 18 energy- 
spectrometer combinations, we find an uncertainty in the 
radii from fits to the entire data set of A te = 0.017 fm, 
ArM = 0.025 fm. These are roughly 2.5 times larger than 
the values quoted in the Mainz analysis, due mainly to 
the increase from their ~0.2% to our 0.5%. Other func¬ 
tional forms give similar results, with the largest effect 
coming from scaling the uncertainties with 1 /Q 2 . The 
cases x = Q 2 , 1 /Q 2 , 9 , and 1/9 are given in Table IX 


We take the case x = 9 to represent a reasonable average 
of the functional forms tested. 


3. Impact of applying systematic corrections to different 
data subsets 

Applying the 0.5% correction over the full kinematic 
range for each spectrometer, rather than over the range 
corresponding to a single beam energy, yielded somewhat 
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TABLE IX: Results for changes in the radii under multipli¬ 
cation (top sign) or div ision (bottom sign) by a linear per¬ 
turbation as in Eq. (35 I for each beam energy/spectrometer 
combination, with x = Q 2 , l/Q 2 , 8, or 1/6. Fits are for the 
657 point rebinned A1 MAMI dataset with 0.3%-0.4% un¬ 
correlated systematic uncertainties at three values of Q 2 lax 
using the z expansion with to = 0, Gaussian priors with 
|Cfc|max — |5fc|max/Pp — 5, fcmax — 12. 


X 

Qmax [GeV 2 ] 

Ar E [fm] 

Ar M [fm] 

Q a 

0.05 

=F0.017 

±0.021 


0.5 

3=0.016 

3=0.022 


1 

3=0.015 

3=0.026 

1 IQ 1 

0.05 

±0.041 

3=0.046 


0.5 

±0.025 

±0.016 


1 

±0.023 

±0.021 

8 

0.05 

3=0.022 

±0.027 


0.5 

3=0.018 

3=0.021 


1 

3=0.017 

3=0.025 

1/8 

0.05 

±0.036 

3=0.039 


0.5 

±0.024 

±0.018 


1 

±0.021 

±0.022 


smaller uncertainties for re and somewhat larger uncer¬ 
tainties for vm■ There was also a wider spread in the 
uncertainties arising from different functional forms in 
Eq. (35), as expected for the interpolation over a wider 
kinematic range. Applying the 0.5% variation only over 
the angular range for each normalization subset yielded 
uncertainties that were typically 20%-30% larger for te 
compared to the default approach, with smaller increases 
for the uncertainty on vm- We note that similar studies 
using the original 1422 point data set showed much larger 
increases when applying the correction to the different 
normalization subsets. 

For simplicity, we have taken the systematic scaling 
factor, a in Eq. ( |35| , identical in sign (i.e., always mul¬ 
tiplying or always dividing by 1 + <5 corr ) and magnitude 
for each data subset. However, many systematic effects 
could differ for the different spectrometers, and the com¬ 
bined effect might be enhanced or suppressed by the as¬ 
sumption of identical corrections. When applied individ¬ 
ually to each spectrometer, the charge radius uncertainty 
tends to be dominated by the corrections applied to spec¬ 
trometer B. For the magnetic radius, there tends to be 
a significant cancellation between the corrections from 
the three spectrometers, and the result of shifting all of 
the spectrometers identically (used in Ref. (5] and shown 
in Table IXI is much smaller than the result of evalu¬ 
ating the corrections independently for each spectrome¬ 
ter. Because it is not clear how much the spectrometer 
corrections may be related, we do not enhance the uncer¬ 
tainty in tm- We simply note that the uncertainty on r« 
shown in Table |IX| could be a significant underestimate 
if the cancellation in these tests does not reflect the true 
nature of any systematic corrections. 

To further investigate the impact of applying differ¬ 
ent correlated systematic shifts to different data subsets, 


consider a fit with distinct parameters a in Eq. (35) for 
different data subsets. These are allowed to vary as part 
of the fit, which then allows for subpercent kinematic 
variations in each data subset. This could be done with 
separate parameters for each spectrometer, for each of 
the 18 energy-spectrometer combinations or for each of 
the 34 different normalization subsets. For definiteness, 
we consider a fit with an independent normalization and 
slope parameter a for each of the 34 normalization sub¬ 
sets. 18 This seems most consistent with the breakdown 
of uncertainties into normalization, correlated systemat- 
ics, and uncorrelated systematics. For Q 2 lax = 0.5 GeV 2 , 
we find te = 0.891(18) fm and tm = 0.792(49) fm, com¬ 
pared to te = 0.895(20) fm and rju = 0.776(38) fm from 
Table |X| below, which includes both statistical and uncor¬ 
related systematic uncertainties. The changes in the ex¬ 
tracted radii are consistent with the previously assigned 
uncertainties associated with the correlated systematics. 
The uncertainties in this fit are somewhat smaller for the 
charge radius and larger for the magnetic radius, in line 
with the expectation based on applying the corrections 
separately to each spectrometer. This may be a more 
realistic estimate of the uncertainties and could poten¬ 
tially allow for a combined analysis of Mainz and world 
data by including all of the Mainz systematic uncertain¬ 
ties explicitly in the fit. However, most likely neither the 
Mainz assumption that the corrections are totally corre¬ 
lated between settings nor the assumption here that they 
are totally uncorrelated is entirely realistic. The analysis 
presented here is included only as an independent esti¬ 
mate of the impact of allowing the correlated systematic 
correction to differ for different kinematic settings. 


/. Final evaluation of the correlated systematics 

It is difficult to determine an optimal approach for 
evaluating the impact of unknown systematic errors or 
corrections. The analysis strategy for the Mainz data 
set involves a breakdown of the uncertainties into uncor¬ 
related, correlated, and normalization contributions and 
seems most consistent with applying the correlated un¬ 
certainty to each normalization subset. As noted above, 
this tends to increase the uncertainty on rE in com¬ 
parison to applying the correlated uncertainty to each 
spectrometer or to each spectrometer-energy combina¬ 
tion. Similarly, applying corrections independently to the 
three spectrometers tends to decrease the uncertainty on 
te and increase the uncertainty on vm- 

We choose to evaluate the correlated systematic error 
by making simple, minimal changes to the A1 procedure. 


18 Note that in the Mainz analysis and our other fits there are 31 
normalization parameters which appear in 34 different combina¬ 
tions. For this test, we allow all 34 normalization factors to vary 
independently. 
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We evaluate the impact of a linear angle-dependent cor¬ 
rection (a; = 0 ), applied to each beam-spectrometer com¬ 
bination, but choose a fixed 0.4% variation. The 0.4% 
variation (a = 0.004) is approximately twice the typical 
value considered in the A1 analysis, which seems a rea¬ 
sonable estimate to account for additional systematic ef¬ 
fects such as TPE [T2] and background subtraction m- 
At Qm ax = 1 GeV 2 , this choice yields uncertainties of 
0.014 fm and 0.020 fm for ve and tm, respectively, four- 
fifths of the uncertainties shown for x = 8 in Table IIXI 
These uncertainties are significant, but not sufficient 
to explain the discrepancy with muonic hydrogen. Ob¬ 
taining larger shifts due to such corrections would require 
either a systematic shift above the 0.4% assumed here, 
a correction applied over smaller data subsets, a more 
extreme functional form for such corrections than con¬ 
sidered here, or a conspiracy between shifts applied to 
different spectrometer-beam combinations. 


VII. RADIUS RESULTS FROM MAINZ AND 
WORLD DATA 



QL* [GeV 2 ] 


Having completed our systematics studies, we proceed 
to perform a final fit to the Mainz data and compare 
with a fit to other world data using the same theo¬ 
retical framework. We close this section with several 
consistency checks on the fits, including a discussion of 
form factor priors, radiative corrections beyond TPE, 
and the verification of the fit consistency between dif¬ 
ferent spectrometer-energy subsets of the data. 


A. Best fit radii from Mainz data 


Let us summarize our final fit to the Mainz data set. 
We use the 657 point rebinned data set of Sec. |VIB 2[ 
with the SIFF sum of monopoles TPE correction from 
Table [I] in place of the Feshbach correction applied in 
the Mainz analysis and with the Al statistical rescal¬ 
ing replaced by a fixed 0.3%-0.4% uncorrelated system¬ 


atic as in Sec. VIB 3 We employ our default form fac¬ 


tor scheme, to = 0, k max = 12, and Gaussian bound 
M max — Mmax/Vp = 5. The results are shown in Fig. [6] 
and Table [X] The “statistical” uncertainty accounts for 
both counting statistics and the uncorrelated systematic 
uncertainties. The energy cut correction is taken from 
Table |VTTT] The correlated systematic uncertainty is ob¬ 
tained from the x = 9 entry of Table IX rescaled to 0.4%, 
as described above in Sec. IVI C 41 


B. Best fit radii from world data 

Now that we have a procedure for analyzing the Mainz 
data, we perform a similar fit to the global set of world 
data, excluding the Mainz data set. We perform this 
separate analysis in part to obtain independent results 


FIG. 6: Extracted electric (top panel) and magnetic (bot¬ 
tom panel) radii as functions of the kinematic cut ax on 
the momentum transfer for the rebinned 657 point Al MAMI 
data set with 0.3%-0.4% uncorrelated systematic uncertain¬ 
ties, using the 2 expansion with to = 0, Gaussian priors with 
| Q’k |max = |6fc|max/p P = 5, fe ma x = 12. Error bands are statis¬ 
tical and uncorrelated systematic only. 


TABLE X: Final radius results from the fits to the rebinned 
657 point Al MAMI data set with 0.3%-0.4% uncorrelated 
systematic uncertainties in Fig. [6] for three values of Q^ ax . 
The uncertainties labeled “stat” include both the statistical 
and uncorrelated systematic uncertainties, while those labeled 
“AE” and “cor” account for the energy cut dependence and 
the correlated systematic uncertainties, respectively. 


Qmax 

rs 

VM 

[GeV 2 ] 

[fm] 

[fm] 

0.05 

0.856(27) s tat(5) A E(18)cor 

1-11(14) st at (2 )a E (2) C or 

0.5 

0.895(14) s tat (3) AE (14)cor 

0.776(34)stat(3)AE(17) cor 

1 

0.908(13)stat (3) AE(14)cor 

0.766(33)stat(9)AE(20)cor 


as a check on consistency between the Mainz data set 
and the world data set. In addition, it is not clear that 
there is a reliable way to perform a combined analysis of 
the Mainz data with other experiments, given the very 
different manner in which uncertainties from the Mainz 
experiment were presented m- The inclusion of the cor¬ 
related systematic correction coefficients a in the fit, as 
discussed in Sec. |VI C 3[ yields a fit to the Mainz data 
where all uncertainties are accounted for in the fit and 
would allow a combined analysis with the world data. 
However, this approach allows the correlated systematic 
corrections to be different for each subset, which may not 
be significantly better than the assumption in the Mainz 
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TABLE XI: Final radius results from the fits to the world 
cross section data in Fig. [7] (first line) and to the combined 
world cross section and polarization data in Fig. [8] (second 
line). There is no polarization data below Q 2 = 0.05GeV 2 . 


G 2 

max 

(GeV 2 ) 

r E 

(fm) 

rM 

(fm) 

x 2 

N a 

N r at 

N gxp 

0.05 

0.846(42) 

1.04(11) 

52.9 

111 

0 

8 

0.5 

0.910(25) 

0.919(38) 

163.4 

269 

0 

15 


0.927(24) 

0.899(38) 

234.5 

269 

30 

15 

1 

0.916(24) 

0.914(34) 

260.9 

363 

0 

23 


0.919(23) 

0.913(34) 

366.0 

363 

41 

23 


analysis that these corrections are identical for different 
subsets. We thus present separate fits to the Mainz and 
world data sets so that the comparison can be made with¬ 
out worrying about how to consistently treat Mainz and 
world uncertainties. 

For the analysis of world data, we take the % 2 function 
Xw = xl + Xb + xl ■ (36) 


Here, y 2 and y 2 are identical to those used for the Mainz 
analysis, Eqs. (32) and (33). Because these experiments 
provide a normalization uncertainty, we follow previous 
analyses and include y 2 for the floating normalization 
parameters assigned to each experiment, where 


X 


2 

n 


f’exp 

E 


(1 - Vi, fit) 2 

dvf 


(37) 


Below Q 2 = 1 GeV 2 , there are N exp = 23 independent 
experiments, and drji ranges from 1.5% to 4.6%. For the 
fit to the combined world and polarization data sets, we 
include an additional term for the recoil polarization and 
polarized target measurements of n p G E / G m , 


,2 _ ,2 
Xw+p x w 


yV {Ri — Ri, fit) 2 

i=l 


dRj 


(38) 


where Ri = n P G E {q 2 )/G M (q?)- 
The TPE model for the cross section data in these 
analyses is the SIFF prescription with the form factor 
as a sum of monopoles from Table | 19 As in the fit to 
the Al MAMI data set, we find little difference in the 
r E results using either a dipole form factor or the Fesh- 
bach correction but significant differences in the 7'm re¬ 
sults between approaches. For Q^ ax = 1 GeV 2 , there 
are N a = 363 cross section data points, Wat = 41 po¬ 
larization data points, and JV exp = 23 normalization pa¬ 
rameters. Results for fits using the 2 expansion with 


19 The use of different conventions in the world data to isolate the 
IR finite TPE contribution, as detailed after Eq. i ] *23 [ l. changes 
r E and r_\,j by less than 0.003 fm. 




FIG. 7: Extracted electric (top panel) and magnetic (bot¬ 
tom panel) radii as functions of the kinematic cut ax 
on momentum transfer for the world cross section data set, 
using the z expansion with to = 0, Gaussian priors with 
|afc|max = |6fe|max/Mp = 5, Wax = 12. Error bands include 
statistical and systematic uncertainties. 


our default to = 0, Wax = 12, and Gaussian bounds 
|afe|max = IMmax/w = 5, are displayed in Figs. [7]andj8j 
Table |XI| contains the radii values and error bud get for 
particular values of Qjhax- The results in Table [Xl] in¬ 
dicate that the inclusion of the polarization data does 
not significantly change the extracted radii. The results 
for the electric radius are in agreement with the fit to 
the Mainz data in Table [Xj while magnetic radius values 
disagree by 2.7er if the uncertainties are added in quadra¬ 
ture. 


C. Consistency checks 


We have derived best-fit values for r E and from 
the Mainz data set and from a world data set excluding 
the Mainz data. We observe a significant dependence of 
the Mainz radius on the Q 2 range included in the fit as 
well as a disagreement between the Mainz and world data 
extractions of tm- Here, we describe several consistency 
checks on the fits we have performed. We also consider 
the possibility of a common systematic not specific to a 
particular experiment, and reexamine subleading radia¬ 
tive corrections which become enhanced at large Q 2 . 
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FIG. 8: Same as Fig. [7] but with both the world polarization 
data in addition to the world cross section data. 


1. Priors 

Let us revisit the dependence on the class of form 
factors over which the fit is performed, defined in the 
bounded z expansion by the choice of k max , to, and coef¬ 
ficient bound. 

As discussed above in Sec. |VD 1[ we have taken & max 
large enough such that the fit results are independent of 
the precise value of /c max , removing this choice from the 
discussion of prior dependence. 

With our imposition of coefficient bounds, the chosen 
form factor class depends on to- 20 We have redone se¬ 
lected fits with different scheme choices, e.g., to = tg pt 
defined after Eq. ([8]) , finding negligible dependence on t 0 
in the large fc max limit. 

Regarding Gaussian vs sharp priors, we have employed 
Gaussian priors for numerical ease but have checked that 
our results are not significantly changed if sharp priors 
are used. Central values for both i'e and r« differ by a 
negligible amount between the two priors, and the differ¬ 
ence in radius errors is small. 

We note that enforcing a bound on the radius param¬ 
eters could in principle bias the radius fits. For example, 
at t 0 = 0 (or default choice), the squared radii are pro¬ 
portional to fit parameters ai and &i, and a bound on 


20 Modifications (multiplication by suitable analytic function (ft) to 
the z expansion can ensure that the form factor class is rigorously 
independent of to in the large fc max limit if the bound is placed 

on i2k a t' See ’ e -§-’ Ref - E 3- 


these parameters would tend to bias fits toward smaller 
radii. We have checked that fitting with the bounds on 
di and bi removed has a negligible impact. 

Finally, consider the choice of the numerical value for 
the bound. We choose a Gaussian bound of 5, i.e., 
|afc/ao|max = |frfc/&o|max = 5, for our fit s, base d on the 
sum rules and studies discussed in Sec. cub Our im¬ 
plementation of the bounds is meant to be very conser¬ 
vative, especially at large fc, where the coefficients must 
fall as 1/fc 4 . We examined the impact of tightening the 
constraint for larger values of fc, taking a bound of 5 for 
k = 1,..., 4 and a bound of 20 /k for larger k values. This 
yields small changes in the central radius values (< 0.3 cr), 
with only slightly smaller uncertainties. Because tighter 
high-A: constraints have a minimal impact, for simplicity, 
we use a fixed bound of |a fc /a 0 | m ax = |&fc/^o|max = 5. 
More aggressive priors or fc max truncations could be in¬ 
voked to reduce the statistical/fit uncertainty on the ra¬ 
dius at the expense of introducing model-dependent trun¬ 
cation errors. This may allow for the possibility of re¬ 
duced uncertainties in the extracted radii, if one can ver¬ 
ify that the reduction in uncertainty coming from tighter 
bounds or truncations is not replaced with a bias that 
yields a larger net uncertainty. In this work, we are fo¬ 
cused on minimizing any such biases and so do not at¬ 
tempt to further constrain the fits. 

We also fit the data and obtained statistical errors 
for larger bounds, |a fc /a 0 | m ax = |^fc/&o|max = 10. The 
change in the extracted radii is very small for fits with 
large & max (in particular for our default fc max = 12). No 
additional uncertainty is applied as typical changes in the 
fit were significantly smaller than the statistical or the 
correlated systematic uncertainties. Fits in which & max 
was not sufficiently large to give fully converged results 
showed larger changes, but are not included in the final 
results presented here. 


2. Data set exclusions 


To verify that the fits to the Mainz data set are not 
biased by one particular subset of the data, we redo our 
best fits for te and Tm 18 times, excluding in each fit 
a particular energy/spectrometer combination. The re¬ 
fer QLx = 0.5 


XIII 


suits are displayed in Tables XII and 
and 1 GeV 2 , respectively. For the electric radius, the im¬ 
pact of each subset exclusion is typically less than half of 
the total statistical error (taken from Table 0- Several 
subset exclusions impact the magnetic radius at a level 
comparable to or greater than the total statistical error. 
For the 180 MeV data, excluding any one of the three 
spectrometers gives a 0.030-0.041 fm shift in r«. This 
is much larger than the estimated systematic uncertainty 
and much larger than one might expect based on an ex¬ 
clusion of 4%-ll% of the data points. While this suggests 
the need for a larger uncertainty in the quoted value of 
rjf, it is hard to quantify what uncertainty would be ap¬ 
propriate as we are comparing highly correlated fits. As 
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TABLE XII: Change in extracted rg and tm when each 


data subset is excluded. Fits are for rebinned A1 MAMI data 
set with 0.3%-0.4% uncorrelated systematic uncertainties at 
Qmax = 0.5GeV 2 (568 data points), using the z expansion 
with t 0 = 0, Gaussian priors with \a k max = |6fc|max//i P = 5, 

k max — 12. 

Spec. 

Beam 

N a 

A r E (fm) 

A r M (fm) 

A 

180 

539 

-0.008 

-0.031 


315 

545 

+0.001 

-0.008 


450 

543 

-0.004 

+0.008 


585 

540 

0.000 

-0.009 


720 

552 

-0.003 

-0.002 


855 

561 

0.000 

0.000 

B 

180 

507 

-0.001 

+0.034 


315 

522 

+0.001 

+0.003 


450 

500 

+0.003 

-0.017 


585 

508 

+0.005 

+0.005 


720 

511 

-0.002 

-0.006 


855 

502 

+0.005 

+0.019 

C 

180 

544 

-0.002 

+0.030 


315 

544 

0.000 

-0.023 


450 

543 

-0.006 

+0.032 


585 

561 

-0.001 

+0.001 


720 

566 

0.000 

0.000 


855 

568 

— 

— 


TABLE 

XIII: Same 

as Table 

XII but for Q^ ax 

= 1 GeV 2 

(657 data points). 




Spec. 

Beam 

N a 

A r E (fm) 

A r M (fm) 

A 

180 

628 

-0.008 

-0.035 


315 

634 

-0.001 

-0.007 


450 

632 

-0.004 

+0.012 


585 

629 

-0.002 

-0.015 


720 

628 

-0.004 

-0.003 


855 

636 

-0.001 

-0.004 

B 

180 

596 

0.000 

+0.041 


315 

611 

0.000 

+0.004 


450 

589 

+0.004 

-0.016 


585 

597 

+0.005 

+0.006 


720 

600 

-0.004 

-0.007 


855 

591 

+0.007 

+0.020 

C 

180 

633 

-0.003 

+0.036 


315 

633 

+0.001 

-0.017 


450 

632 

-0.006 

+0.021 


585 

639 

+0.001 

-0.000 


720 

625 

-0.002 

-0.005 


855 

636 

+0.001 

+0.001 


such, we simply note this as another potential issue, sim¬ 
ilar to the anomalous ax dependence observed in the 
extraction of the charge radius. 


3. Subleading radiative corrections 


We have used standard prescriptions for the electron 
vertex and bremsstrahlung radiative corrections. As 


noted above, in the Q 2 ~ GeV 2 regime, it is critical to 
resum the leading cr log 2 (Q 2 /to 2 ) terms. However, nu¬ 
merically enhanced subleading logarithms can also have 
a significant impact, as illustrated by the following con¬ 
siderations. This exercise is presented both as an illustra¬ 
tion of how a correction would need to deviate from the 
assumptions of Sec. El C 4| in order to reconcile muonic 
hydrogen with electron scattering measurements of the 
charge radius, and to point out the potential impact of 
a class of naively subleading but numerically enhanced 
radiative corrections. 

Recall the explicit form for the sum of the one-loop 
electron vertex and real bremsstrahlung radiative correc¬ 
tions, 


7 T 


, q 2 ; 

log 2 1 

mi 


log 


(r)AE) 2 

EE' 


13 Q 2 
yriog — 
() mi 


(39) 


where r] = E/E' and the ellipsis denotes terms not con¬ 
taining large logarithms. In the regime where A E ~ m e 
(and E ~ E' ~ Q) the numerically relevant leading log 
term 


Q 2 


s - log — 

7 r l mi 


(40) 


is fixed by infrared divergences whose forms are dictated 
by soft photon theorems m- Equivalently, an effective 
theory renormalization analysis between hard (~ Q) and 
soft (~ m e ) scales determines the relevant Sudakov form 
factor. However, in practice, A E can be large compared 
to m e , introducing another scale into the problem, and 
associated large logarithms not captured by the naive ex¬ 
ponentiation of one-loop corrections. A complete analysis 
is outside the scope of the present paper, but to illustrate 
the potential impact, let us consider in place of the ansatz 


that makes the replacement (311 in Eq. (29) the following 
expressions: 


(l + <5) 


l±(d+-log 2 ^ 
7 r mi 


i ±i 


, «, 2 Q 

x exp --log — 
7 r mi 


• (41) 


These expressions agree with the known corrections 
through one-loop order and resum the leading logarithms 
to all orders in perturbation theory when there is only one 
large ratio of scales. 

Figure [9] illustrates the impact of applying the cor¬ 
rection on the right-hand side of Eq. ( 4T|) in place of 
the ansatz ©. For definiteness, the plot takes A E = 
lOMeV. As indicated in the figure, the shifts in the 
radii under this correction are a factor ~ 2-3 larger than 
those allowed in Table m which considered corrections 
varying by 0.5% over beam-energy/spectrometer combi¬ 
nations. The variation of the correction (41) over beam- 
energy/spectrometer combinations [i.e., the magnitude 
of a in Eq. (351] ranges between 0.9% and 2.6%, with an 
average 1.5%% 
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Qmax [GeV 2 ] 


FIG. 9: Illus trat ive fit with modified radiative corrections 
given by Eq. (411 using A E = lOMeV. Lower and upper 
dashed blue lines correspond to the plus sign and minus sign 
in Eq. ( |41| ), respectively. Fits are for the 657 point rebinned 
A1 MAMI data set with 0.3%-0.4% uncorrelated systematic 
uncertainties using the 2 expansion with to = 0, Gaussian 
priors with |a fc | max = |fefc| m ax/Atp = 5, fcmax = 12. Black solid 
lines reproduce the curves in Fig. [6] For orientation, the dash- 
dotted red line indicates the muonic hydrogen value for r_g. 


D. Final radius extractions 


A global analysis combining Mainz and other world 
data will artificially favor the Mainz data, as the uncer¬ 
tainties associated with each cross section measurement 
include only a small part of the total uncertainty. Thus, 
we provide best-fit values separately for our analyses of 
the Mainz and world data. To determine an optimal 
Qmaxi Fig. 10 illustrates the statistical uncertainty on 
te and tm found using our default fit both to the 1422 
point Mainz data set and to the world data set. For 
the Mainz data, the uncertainty is minimized by taking 
Qmax ~ 0.5GeV 2 , with negligible improvement beyond 
this point. To maximize the statistical power of the data, 
while minimizing potential systematic effects in higher 
Q 2 data, we take for definiteness the Qmax = 0.5 GeV 2 
results of the previous sections. 21 


21 A similar choice was made in Ref. m based on radius sensitivity 
in the world data summarized by extracted form factors [26]. Re¬ 
lated argumentation for the ax dependence of radius sensitiv¬ 
ity, based on continued-fraction expansion, is given in Ref. m- 



FIG. 10: Statistical error on te (bottom, red squares) and tm 
(top, blue circles) as a function of Q alax . Solid symbols are 
for the 1422 point A1 MAMI data set, and open symbols are 
for the world cross section and polarization data set. Fits use 
the z expansion with to = 0, Gaussian priors with |ofe| m ax = 
|5fc|max/ /Tp — 5, femax — 12. 


We then have for the Mainz data set, from Table [Xj 


Mainz = o.895(14)(14) , r^ ainz = 0.776(34)(17) , (42) 


' E 


where the first error comes from counting statistics and 
uncorrelated systematics and the second error comes 
from variation of the bremsstrahlung energy cut and cor¬ 
related systematics. For the world data set, including 
cross section and polarization measurements, w e ta ke a 
slightly higher Q;^ ax = 0.6 GeV 2 based on Fig. 10 We 
then have 


T world = 0 _ 916 ( 24 ) , 


world 

r M 


= 0.914(35). 


(43) 


These values c orre spond to the same analysis as pre¬ 


sented in Table 


XI 


but for a fit with Q, nax = 0.6 GeV . 


In contrast to the Mainz data, the world data have com¬ 
bined statistical and systematic uncertainties at the cross 
section level and so have only a single combined uncer¬ 
tainty. 

The electric charge radius results are consistent with 
each other and between one and two standard devia¬ 
tions higher than the atomic physics measurements based 
on atomic hydrogen which yield te = 0.8758(77) fm [I]. 
They are well above the muonic hydrogen result 7'e = 
0.84087(39) fm [4]. The magnetic radii differ significantly, 
indicating an unresolved tension between the Mainz data 
set and the world data set. 


A simple combination of the results (421 and (43) yields 


r^ vg - = 0.904(15), r^ g - = 0.851(26). (44) 
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While the Mainz and world data sets have comparable to¬ 
tal uncertainties, the high statistics of the Mainz dataset 
imply that in this case the errors are dominantly sys¬ 
tematic. It is not entirely clear that a simple average of 
the Mainz and world results is appropriate mm - The 
magnetic radii differ by 2.7 ct, suggesting an inconsistency 
between Mainz and world cross sections that is ignored 
in averaging the results. In addition, the simple combi¬ 
nation assumes that the uncertainties for the Mainz and 
world data analyses are independent, which may not be 
the case if there is a common error, e.g., due to approxi¬ 
mations in the radiative correction procedures. 


VIII. SUMMARY AND DISCUSSION 

We have performed a comprehensive analysis of 
electron-proton scattering data to determine the proton 
electric and magnetic radii. Our analysis incorporates 
constraints of analyticity and perturbative scaling which 
enforce model-independent bounds on form factor shape. 
The bounded z expansion ensures that the true form fac¬ 
tor is guaranteed to lie within the space of considered 
curves, while at the same time being sufficiently restric¬ 
tive to enable meaningful radius extractions. We focused 
on the high-statistics Mainz data set, and performed a 
wide- ranging study of the impact of potential systematic 
errors. We discussed potential flaws in the procedure of 
rescaling statistical errors and addressed these by rebin¬ 
ning data taken at identical kinematic settings and apply¬ 
ing a constant uncorrelated systematic error that is not 
assumed to scale with statistics. We also reevaluated the 
correlated systematic uncertainties, increasing the size 
of these effects to include contributions neglected in the 
original analysis, and examining different approaches to 
evaluating the impact of such corrections on the radius. 

Table |XIV| displays the progression of results leading 
up to the final Mainz radius values, as various improve¬ 
ments are included in the analysis. The data exhibit 
several issues that suggest the need for additional uncer¬ 
tainties, but for which it is difficult to quantify an ap¬ 
propriate correction or uncertainty contribution. There 
is an unusually large variation of te and tm with the Q 2 
range included in the fit, as illustrated in Fig. [6] In addi¬ 
tion, the exclusion of individual data sets, in particular 
at low beam energy, has an unusually large impact on 
the extracted radii. Despite these anomalous features, 
inclusion of the above improvements leads to a proton 
charge radius that is larger than extracted in the original 
A1 analysis |9], and, even with the larger uncertainty, 
almost 3er above the value r e ~ 0.84 fm inferred from 
muonic hydrogen. Based on our examination of the sys¬ 
tematic uncertainties, resolving this discrepancy would 
require correlated systematic effects well above the 0.5% 
level that was considered in our analysis. 

As an independent check of the radius, we performed 
the analogous fit to the world data excluding Mainz data. 
The systematic error treatment in the world data set dif- 


TABLE XIV: Charge and magnetic radii as determined in 
Ref. [9] compared to the sequence of fits leading to the fi¬ 
nal values determined in this paper. For the Mainz data set, 
the first error is a combination of statistics and uncorrelated 
systematics, and the second error is from correlated system- 
atics. The entry labeled “alternate approach” is the test from 
Sec. El C 3 1 which evaluates the impact of the correlated sys¬ 
tematic uncertainties as part of the fit, rather than evaluating 
it separately. 


Source te (fm) 


I'm (fm) 


A1 spline 
Bounded z exp. 
+Hadronic TPE 
Rebin, 0.3%-0.4% syst. 
+0.4% corr. syst. 

Ql ax = 0.5 GeV 2 
(Alternate approach) 


Tab 

Tab 


L 


III 


Tab 
Tab IX 


(421 


vr 


0.879(5)(10) 

0.920(9%-) 

0.918(9%-) 

0.908(13%-) 

0.908(13%14) 

0.895(14%14) 

0.891(18%-) 


0.777(13)(14) 
0.743(25) (-) 
0.780(25) (-) 
0.767(33) (-) 
0.767(33) (22) 
0.776(34) (17) 
0.792(49) (-) 


New fit to world data ( 

43 

1 0.916(24) 0.914(35) 

Simple avg. ( 

44 

1 0.904(15) 0.851(26) 


fers in that the systematic errors are included in each 
cross section data point as opposed to deduced from a 
combination of statistics rescaling and model-dependent 
correlated systematics analysis. It is thus not straightfor¬ 
ward to perform a meaningful combined fit, but observ¬ 
ables such as radius values may be compared to verify 
consistency. In this comparison, the Mainz and world 
rE values are in good agreement, and the rn values dif¬ 
fer by 2.7cr. This is perhaps not surprising, given the 
clear disagreement between the Mainz form factors and 
world data, in particular for Gm at low Q 2 [S], Putting 
aside the discrepancy in magnetic radii, the charge ra¬ 
dius puzzle persists, with the world value for r# 3cr high 
compared to muonic hydrogen and the combined Mainz 
+ world average discrepant at the 4.2cr level. In light of 
this, it is also important to inquire to what extent an un¬ 
derestimated systematic effect or theoretical correction 
could be common to both data sets. 

The theoretical input to the radius extraction consists 
of specifying the form factor class and defining the radia¬ 
tive correction model. We have examined in detail the 
uncertainties associated with form factor shape assump¬ 
tions. We find a large impact from fitting to the phys¬ 
ical form factor class defined by the bounded z expan¬ 
sion, compared to polynomial or inverse polynomial fits. 
Somewhat surprisingly, the central value for the charge 
radius goes in the direction of increasing tension between 
the electron scattering and muonic hydrogen. We have 
further examined the dependence of radius values on form 
factor priors, finding that such a residual dependence is 
small compared to other uncertainties. 

The other theoretical input is the radiative correction 
model, as described in Sec. |IVj For the most part, the 
corrections are known precisely or, for model-dependent 
terms including hadronic vacuum polarization or proton 
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vertex corrections, the uncertainties are estimated to be 
small compared to the uncertainty in the radius extrac¬ 
tion. Through one-loop order, the only essential model 
dependence of the radiative corrections enters from the 
description of the TPE. Varying over models in the lit¬ 
erature reveals no large dependence on the applied TPE 
correction. The Q 1 2 3 ~ GeV 2 regime demands QED radia¬ 
tive corrections beyond one-loop order. In the counting 
m e ~ A E, leading logarithms are resummed by a stan¬ 
dard ansatz. Subleading logarithms then enter at a level 
expected to be contained within the 0.4% systematic er¬ 
ror budget. Possible enhancements, either simply numer¬ 
ical or due to the hierarchy A E S> m e , could potentially 
give rise to larger effects. The Mainz and world data 
sets differ in their treatment of bremsstrahlung radiation 
and approximations based on (311, and in the uncertain¬ 
ties ascribed to these effects. More refined calculations, 
including a detailed examination of experimental condi¬ 
tions and the interplay with background modeling and 
subtraction, are required in order to fully address this 
question m 

Further constraints may be placed on the proton form 
factors in combination with electron-proton scattering 
data. In particular, the inclusion of either electron- 
neutron scattering data, or both electron-neutron and 
pion-nucleon data, allows the threshold t cut appear¬ 
ing in the definition of the z expansion (Jsj) for the 
isoscalar/isovector decomposition of the form factors to 
be raised from 4m 2 to either 9m 2 or 16 m 2 . This yields a 
smaller |z| max and hence tighter constraints on the form 
factors and smaller radius uncertainties. Isospin violating 
corrections and systematic uncertainties in the additional 
data must be properly accounted for. These additional 
constraints by themselves cannot offer a satisfactory reso¬ 
lution to the proton radius puzzle, since they would then 
be inconsistent with the results of the fit to the electron 
scattering data alone. Similar remarks hold for the model 
spectral function analysis in Ref. m ■ We note that it 
is not feasible to reconstruct accurate spectral functions 
for form factors, ImG(f), from scattering data, since these 
have support at |z| = 1. 

Electron scattering data from a polarized target have 
been taken and will provide low-Q 2 measurements of 
H p Ge/Gm, down to Q 2 « 0.02 GeV 2 , yielding a more 
sensitive measurement of the magnetic form factor in 
this region mi Future experiments will obtain low-Q 2 
(~ 10~ 4 5 — 10 -2 GeV 2 ) proton form factor measurements 


less prone to high-Q 2 systematics fSOJIHI]- Muon-proton 
scattering promises to provide further insight [82]. In¬ 
dependent of scattering measurements, new results are 
anticipated from hydrogen spectroscopy that may impact 
the charge radius discrepancy, including a new microwave 
measurement of the 2S-2P Lamb shift [83], measure¬ 
ment of 2S-4P transitions [M], and 1S-3S/3D transi¬ 
tions [551 186] . Next-generation lattice QCD simulations 
may provide another handle |87l - 191j : in particular, re¬ 
solving the ~ 8% discrepancy between te ~ 0.84 fm of 
muonic hydrogen and ss 0.91 fm of the simple fit to 
Mainz electron scattering data is a viable present-day 
target. Such first-principles calculations would be inde¬ 
pendent of radiative corrections involving the electron, 
thus avoiding the reliance on hadronic models for TPE, 
and detector-dependent modeling of radiative tails. 

Regardless of its resolution, the proton radius puzzle 
has important implications across particle, nuclear, and 
atomic physics. For example, understanding and control¬ 
ling systematic effects, including radiative corrections, at 
the percent level will be critical for measurements at fu¬ 
ture long baseline neutrino experiments [921 . Any de¬ 
ficiency in the theoretical treatment of electron-proton 
scattering will be exacerbated in neutrino applications 
by the presence of additional flux uncertainties and nu¬ 
clear corrections. Much more intriguing is the possibil¬ 
ity that updated measurements and a detailed exami¬ 
nation of the radiative corrections will not resolve the 
discrepancy. Already much work has been performed to 
find explanations in terms of physics beyond the Stan¬ 
dard Model [7[ |5] I55H55] . Future measurements will pro¬ 
vide more stringent tests of the discrepancy in electron 
scattering and atomic hydrogen, with plans to directly 
compare electron-proton and muon-proton scattering as 
a test of lepton nonuniversality. 
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